- Email: [email protected]
Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank
YMGME-05929; No. of pages: 6; 4C: Molecular Genetics and Metabolism xxx (2015) xxx–xxx
Contents lists available at ScienceDirect
Molecular Genetics and Metabolism journal homepage: www.elsevier.com/locate/ymgme
Regular Article
Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank Jonas Grauholm a, Sok Kean Khoo b, Radoslav Z. Nickolov c, Jesper B. Poulsen a, Marie Bækvad-Hansen a, Christine S. Hansen a, David M. Hougaard d, Mads V. Hollegaard a,⁎ a
Section of Neonatal Genetics, Danish Centre for Neonatal Screening, Department of Congenital Diseases, Statens Serum Institut, Artillerivej 5, DK-2300 Copenhagen S, Denmark Department of Cell and Molecular Biology, Grand Valley State University, Grand Rapids, MI 49503, USA Department of Mathematics & Computer Science, Fayetteville State University, Fayetteville, NC 28301, USA d Danish Centre for Neonatal Screening, Department of Congenital Diseases, Statens Serum Institut, Artillerivej 5, DK-2300 Copenhagen S, Denmark b c
a r t i c l e
i n f o
Article history: Received 5 May 2015 Received in revised form 26 June 2015 Accepted 26 June 2015 Available online xxxx Keywords: Biobank Dried blood spots Neonatal screening RNA storage Gene expression profile Whole transcriptome amplification
a b s t r a c t A large part of the human genome is transcribed into various forms of RNA, and the global gene expression profile (GEP) has been studied for several years using technology such as RNA-microarrays. In this study, we evaluate whether neonatal dried blood spot (DBS) samples stored in the Danish Neonatal Screening Biobank (DNSB) can be used for GEP. This paper is divided into sub-studies examining the effects of: 1) different whole transcriptome amplification kits (WTA); 2) years of storage and storage in room temperature (RT) versus freezers (−20 °C) on DNSB DBS samples; 3) effects of RT storage vs freezer storage on DBS samples from the USA and DNSB, and 4) using smaller disc sizes, thereby decreasing DBS use. We present evidence that reliable and reproducible GEPs can be obtained using neonatal DBS samples. The main source of variation is the storage condition. When samples are stored at −20 °C, the dynamic range is increased, and Pearson correlations are higher. Differential analysis reveals no statistically significant differences between samples collected a decade apart and stored at −20 °C. However, samples stored at RT show differential expression for a third of the gene-specific probes. Our data also suggests that using alternate WTA kits significantly changes the GEP. Finally, the amount of input material, i.e., the size and number of DBS discs used, can be reduced to preserve this valuable and limited material. We conclude that DNSB DBS samples provide a reproducible resource for GEP. Results are improved if the cards are stored at −20 °C. Furthermore, it is important to use a single type of kit for analysis because using alternate kits introduces differential expression. © 2015 Elsevier Inc. All rights reserved.
1. Introduction RNA is a diverse molecule essential to cellular function. One form carries the genetic code from the nucleus while others control diverse regulatory functions, such as splicing and translation. The ENCODE project showed that up to 93% of the human genome can be transcribed into one of several known forms of RNA [1–3]. Since the introduction of genome wide association studies, approximately 2000 studies have linked numerous diseases to more than 10,000 SNPs. The challenge now is to understand their functional consequences. Especially in complex diseases, it is hard to infer organism level ⁎ Corresponding author at: Section of Neonatal Genetics, Danish Center for Neonatal Screening, Department of Congenital Diseases, Statens Serum Institut, Artillerivej 5, DK-2300 Copenhagen S, Denmark. E-mail addresses: [email protected] (J. Grauholm), [email protected] (S.K. Khoo), [email protected] (R.Z. Nickolov), [email protected] (J.B. Poulsen), [email protected] (M. Bækvad-Hansen), [email protected] (C.S. Hansen), [email protected] (D.M. Hougaard), [email protected] (M.V. Hollegaard).
phenotypes (i.e., disease state) from observed genotypes. Changes in molecular phenotype, as observed through up- or down-regulation of RNA, may help bridge the gap between genotypes and phenotypes [4–6]. Numerous technologies and methodologies are used to study RNA, and studies of gene expression profiles (GEP) are expected to have an increasing impact in research and as a clinical tool [7–9]. Denmark has a long tradition of collecting and storing data in nationwide registries. The registries, which can be linked through a unique person identifier number (CPR-number), contain information from birth to death, including notable social and health-related events [10]. The accessibility of well-defined samples is generally a limiting factor in genetics. In addition to the data registries, Denmark also has a long history of storing biological material, including the excess neonatal dried blood spot (DBS) samples routinely collected as a part of the Danish Neonatal Screening program. The DBS samples, stored in the Danish Neonatal Screening Biobank (DNSB), are primarily used to improve neonatal screening and diagnostics [11]. It is also possible to access the DBS samples for research on improving the Danish population's
http://dx.doi.org/10.1016/j.ymgme.2015.06.011 1096-7192/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
2
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
health in general. By January 2015, the DNSB contained more than 2.2 million DBS samples dating back to 1982. Combining the biological material from the DNSB DBS samples with data from the registries unlocks a powerful tool for conducting unbiased population-based health research which are based on well-documented diagnoses and thorough anamnesis. We have previously shown that DNSB DBS samples can be used for array genotyping [12,13], DNA methylation [14], and whole genome and exome sequencing [15]. Using these approaches, the DNSB has been successfully used to investigate several diseases [16–21]. Other groups have published RNA work using DBS cards as a source [22–24]. A series of work by Dr. Khoo and colleagues has convincingly shown that DBS samples from the Michigan Screening Biobank (MSB) can be used for GEP [23] and to detect the differential expression of gender specific genes [25]. Dr. Khoo and her colleagues have also documented differences between DBS stored at room temperature versus conventionally preserved blood after more than 10 years in storage [26]. This study evaluates whether DBS samples from the DNSB can be used for reliable GEPs. Moreover, we evaluate the effect of storage conditions, including time and temperature, the effect of using different whole-transcriptome amplification (WTA) kits and whether the GEP remains accurate and stable when using smaller disc sizes, which provide less whole blood for the RNA extraction.
2. Materials and methods
Denmark. The Institutional Review Board at the Michigan State University approved the use of neonatal DBS samples from the MSB. This study consists of four sub-units (see Fig. 1) evaluating the effect on GEPs when: 1) Using two different WTA kits to compare samples collected in the year 2000 and stored at −20 °C. 2) Testing DBS samples stored in different conditions, using two sets of DNSB samples collected in 1982 and 1990. While the 1990 samples were stored at −20 °C shortly after the neonatal screening analysis, the 1982 samples were stored at room temperature (RT) for up to a month before being transferred to −20 °C. Because samples collected in 2000 were stored in the same conditions as the 1990 samples, we also compare the results of this analysis to data generated in Sub-study 1 to determine whether the duration of storage causes variance. 3) Comparing DBS samples stored at RT to those stored at −20 °C. To provide further support that storage conditions, rather than years of storage, are the primary source of variance, we compare technical replicates of neonatal DBS samples collected in 1984 and 2011 and stored at RT in the MSB. As a reference, we use two DNSB neonatal DBS samples collected in 2000 and stored exclusively at −20 °C. 4) Decreasing the amount of whole-blood (WB) used for RNA extraction and, thereby, the WTA RNA input. We analyzed DNSB DBS samples stored at − 20 °C since the year 2000. Three conditions were tested: 2 × 3.2 mm discs (~ 6 μL WB), 1 × 3.2 mm discs (~3 μL WB), and 2 × 1.6 mm discs (~1.5 μL WB).
2.1. Sample material and study design 2.2. RNA isolation, WTA, and microarray technology In this study, we used DBS samples stripped of all information except the year of storage. As such, this study does not constitute a healthrelated research project as defined by the Danish “Act on Research Ethics Review of Health Research Projects”, but is considered a developmental project for the newborn screening program. Development projects for neonatal screening do not require prior approval from the Committees on Biomedical Research Ethics for the Capital Region of
Total RNA was extracted, purified, and concentrated using an Illustra RNAspin mini kit (GE Healthcare, Little Chalfont, Buckinghamshire, United Kingdom) and treated with DNase to eliminate DNA contamination. We used the WT-Ovation Pico RNA amplification system (NuGEN) (NuGEN, San Carlos, California, USA) or the Sigma-Aldrich's complete whole transcriptome amplification kit (Sigma) (Sigma, St. Louis,
Fig. 1. Schematic overview of storage conditions and sample processing of the four sub-units in this study. The MSB neonatal DBS samples were stored at RT prior to the experiment (2 and 30 years respectively). The DNSB neonatal DBS samples are stored predominantly at −20 °C since 1982, 1990 and 2000 (12, 20 and 30 years, respectively). The DBS samples were excised in sizes of 2 × 3.2 mm, 1 × 3.2 mm or 2 × 1.6 mm, their RNA extracted and WTA using NuGEN or Sigma kits, before subjection of the samples to microarray analysis with the Agilent 8 × 60K RNA array. Sub-study 1 compares the effects of WTA kits. Sub-study 2 compares the storage conditions used in the DNSB. Sub-study 3 compare DNSB samples stored at frost with MSB samples stored at RT. Sub-study 4 explores if it is possible to reduce the amount of input material.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
3
Missouri, USA) to generate single-stranded cDNA that was subsequently labeled with the fluorescent dye Alexa Fluor 3 before hybridization onto the Agilent whole human genome gene expression 8 × 60K microarray (Agilent 60K) (Agilent, Santa Clara, California, USA). Each array contains 62,976 oligonucleotide probes covering 27,958 Entrez Gene RNAs and 7419 long intergenic non-coding RNAs. The arrays were hybridized for 17 h at 65 °C and a rotation speed of 10 rpm, washed for 2 min each with wash buffers 1 and 2 (Agilent, Santa Clara, California, USA), and scanned using an Agilent G3 high-resolution scanner (Agilent, Santa Clara, California, USA). Probe features were extracted from the microarray data using Agilent's Feature Extraction software v.10.7.3.1 (Agilent), and used for GEP analysis via Partek Genomics Suite (PGS) version 6.6 (Partek Incorporated, St. Louis, Missouri, USA). 2.3. Data analysis The gProcessedSignal and gIsInNegCtrlRange values were loaded into PGS. The gProcessedSignal measurement is an analogue value based on a given probe's intensity and is assumed equivalent to the expression level of the transcript corresponding to the probe. The gIsInNegCtrlRange measurement is a digital value that is set as “1” when the intensity value of the gProcessedSignal is within three standard deviations of the array background. Probes that were not flagged (i.e., those with gIsInNegCtrlRange values of 0) are considered detected. The Initial QC was based on the gProcessedSignal value for the unfiltered set of ~60K probes. Evaluation through principal component analysis (PCA) was based on a covariance model and hierarchical clustering analysis (HCA) based on a Euclidean model of complete linkage. To generate Venn diagrams, we quantified the shared and unique entries in the lists of detected probes. We then log2-transformed and quantile-normalized the gProcessedSignal, which was subsequently used to generate intensity–intensity scatterplots of sample correlations and to perform an analysis of variance (ANOVA). The ANOVAs were run as separate one-way analyses, where the factors were the age of the card or the WTA kit. Significant hits were those with at least a 2-fold expression change and a p-value of less than 0.01 following a Benjamini and Hochberg correction for false discovery rate [27,28]. 3. Results and discussion 3.1. RNA microarray quality control The first step in our RNA array quality control analysis was to consult the un-normalized log2 transformed intensities to identify potential biases and samples to exclude. The raw probe intensity QC showed an increased dynamic range for gene-specific probes with higher average intensities compared with the negative control probes (Supplemental File 1); demonstrating that gene-specific probe signals are above background levels in all samples. Not only does this confirm the presence of reliable material in all samples used in the microarray analysis, but it also validates the technical aspects of the experimental setup. Following quantile normalization of the intensity values, we used PGS to perform a PCA using a covariance model. Using the normalized data in PGS, we created an HCA using a Euclidean model of complete linkage. For both methods, we used the complete set of gene-specific probes to avoid selection bias. The PCA and HCA (Fig. 2) show distinctive clustering of the groups. Based on the PCA (Fig. 2A), we were able to define the two main contributors of variance. The storage condition (i.e., − 20 °C versus RT) of DBS samples accounts for 51.9% of the variation (PC1) and the use of different WTA kits accounts for 13% of the variation (PC2). Combined, the two principal components separate samples into groups of four quadrates with non-overlapping separation. The quadrates, separated by the midpoint of the axes, only contain samples from one of the four possible groupings based on the analysis with two WTA kits combined and two storage conditions.
Fig. 2. A) PCA generated using a covariance model. B) HCA constructed using a Euclidian model for complete linkage. Both (A & B) are generated using all gene-specific probes as to avoid selection bias. Each dot in the PCA plot and each line in the HC represent one sample, individual sample names have been removed, instead we have color coded by storage and preparation kit: Sigma-(−20) (green), Sigma-RT (purple), NuGEN-(−20) (red) and NuGEN-RT (blue). Principal component 1 (PC1) stratifies the samples on storage conditions (−20 °C vs RT), whereas PC2 splits by the kit used on the WTA (Sigma vs NuGEN). The HCA shows a similar dependency by WTA kit and storage conditions as in the PCA plot, though the samples stored at room temperature (RT) group as a single branch in the HCA.
Similar to the PCA, the HCA in Fig. 2B shows sample dependency by WTA kit and storage condition. However, the RT storage samples cluster (Fig. 2B). Importantly, when repeating the analysis using only Agilent control probes, neither the PCA nor the HCA produces the pattern shown in Fig. 2 (see Supplemental File 2). This pattern difference suggests that the effect seen in Fig. 2 is the result of variance from storage conditions and the WTA kit used, rather than artifacts introduced by the experimental setup. Our observations are consistent with those previously reported. Storage conditions have been established as a major factor in RNA preservation e.g., in the extreme case, comparing formalin-fixed, paraffinembedded samples to their fresh frozen counterparts [29]. Similarly, variances from the use of different WTA kits are well described [30]. In contrast, the duration of DBS sample storage did not play an important role in the observed variances. In both the PCA and HCA, the samples drawn in 1990 cluster with the samples drawn in 2000. Thus, we conclude that storage condition, rather than storage duration, is the primary source of high variance in DBS sample analysis. As a result, we are confident that our measurements reflect the biology of samples rather than an introduced bias; we therefore analyzed the full data set. With respect to storage condition, differences appear between groups, even at the raw data level. In Supplemental File 1, the signalto-noise ratio increases in samples from − 20 °C versus RT storage. Extracting the median signal for all RT or − 20 °C negative controls, we found a probe intensity just above 2 (in both), while the gene-specific probes had average median signals of 3.7 in the RT and 4.8 in the −20 °C, confirming that the signal to noise ratio is increased by storage at −20 °C.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
4
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
3.2. Sub-study 1) whole transcriptome amplification kits Eight samples were analyzed as part of the first sub-study. These were all DNSB samples from the year 2000 that were stored exclusively at −20 °C. Transcription profiles were generated for each sample based on the use of two WTA kits. Consequently, identical samples were processed twice, once using the NuGEN kit and once using the Sigma kit. For each sample, we quantified gene-specific probes unique to one of the kits and those detected in both. The counts were averaged between samples to create a single Venn diagram for all comparisons. We also calculated the Pearson correlations, listed with the Venn diagram in Fig. 3. Fig. 3 shows the cross-kit replications comparing identical samples. The Pearson correlations for all combinations of WTA kit and samples are in Supplemental File 3. Overall, the Pearson correlation is higher between different samples tested with the same WTA kit (Pearson average correlation ± standard deviation for NuGEN: R2 = 0.95 ± 0.00 and Sigma: R2 = 0.97 ± 0.00) than between identical samples amplified using different WTA kits (R2 = 0.91 ± 0.01). The raw counts of unique and commonly detected genes used to calculate the mean values in the Venn diagram (Fig. 3) are shown in Supplemental File 4. We consistently detected a higher number of genes using the Sigma kit versus the NuGEN kit. Next, we did an ANOVA to compare identical sample sets prepared using different WTA kits. After having applied a Benjamini and Hochberg correction [27,28] for false discovery rate, we identified 8462 genes with at least a 2-fold change in expression level and a significance of P b 0.01. These data provide further support that, for experimental groups to remain comparable with respect to specific probe detection and expression levels, it is essential to prepare samples using the same WTA kit.
3.3. Sub-study 2) effects of storage time and temperature The second sub-study further investigates the reproducibility of technical replicates and explores the effect of storage at −20 °C versus RT. This analysis used two sets of DNSB samples: the first collected in 1982 and the second in 1990. Samples collected in 1982 had a latency of up to a month before being transferred to − 20 °C, whereas those collected in 1990 were transferred to −20 °C shortly after screening. We quantified the gene-specific probes detected as unique to one or common to both of the replicates. These were averaged to create a single Venn diagram, which is shown in Fig. 4 along with the Pearson correlations. The Pearson correlations comparing all samples in the sub-study are available in Supplemental File 5. The raw counts of the unique and commonly detected genes used to calculate the means shown in the Venn diagrams in Fig. 4 are in Supplemental File 6.
Fig. 3. Pearson correlations and an averaged Venn diagram for the samples used to compare the WTA kits.
Of the two groups with higher Pearson correlation of technical replicates, there is a significant difference between the 1990 samples (R2 = 0.97 ± 0.01) and the 1982 samples (R2 = 0.91 ± 0.02). Furthermore, Supplemental File 5 shows that the Pearson correlation between any two random samples in the 1990 cohort (R2 = 0.97 ± 0.01) is better than the correlation between technical replicates from the 1982 sample mentioned above. Although we did not study technical replicates for samples from the year 2000, we again observed improved Pearson correlations between unrelated individual samples stored at − 20 °C. In Supplemental File 3, R2 = 0.97 ± 0.00 when comparing any two unrelated samples prepared using the Sigma kit. These data indicate that performance is more consistent for samples directly transferred to storage − 20 °C, i.e., the 1990 and 2000 samples, despite the 10year difference in storage duration. Thus, samples whose collection is separated by decades can be included in the same study as long as storage and experimental conditions are stringently controlled. There is however, a marked reduction in reproducibility when analyzing the 1982 samples, which had been stored at room temperature for up to a month before being transferred to storage at −20 °C for three decades. Pearson correlations for the 1982 sample group are in accordance with previous reports on samples stored at room temperature for longer periods. For similar samples tested on the Agilent 44 k array, Khoo et al. found a correlation for experimental replicates of R2 = 0.92 [23]. Overall, the values reported here, R2 = 0.97 ± 0.01 and R2 = 0.91 ± 0.02 for different sample groups, are above those of experimental replicates reported as satisfactory for FFPE samples (R2 = 0.9) [29]. To further investigate the effect of storage time, we performed an ANOVA analysis on DNSB samples from the years 1982, 1990 and 2000. We did not identify any genes with differential expression while comparing DNSB samples from 1990 with those from 2000. In contrast, we detected 22,356 differentially expressed gene-specific probes of the 58,782 probes tested when comparing samples from 1982 to 1990. A comparison of samples from 1982 and those from 2000 detected 27,660 differentially expressed gene-specific probes. For all comparisons, hits were only retained if the p-value was less than 0.01 following Benjamini and Hochberg correction and if the hit had at least a two-fold change in signal intensity between samples (i.e., a log2 change of ±1). These experiments demonstrate that the RNA in DBS samples stored at −20 °C serves is better preserved, leading to increased success and accuracy in potential RNA array studies. 3.4. Sub-study 3) storage at RT versus −20 °C To further probe the effect of room temperature storage on the reliability of DBS samples used in RNA profiling, we carried out experiments using two MSB samples that had been stored exclusively at room temperatures. To highlight the differences between RT and −20 °C storage, we compared the MSB samples to two DNSB samples. Samples were prepared using the Sigma and the NuGEN kit and, while we have already shown reduced decreased performance in cross kit comparisons from the PCA (see Fig. 2), we expect the storage conditions to be a bigger cause of variation. The Pearson correlations for the MSB samples (RT storage) were R2 = 0.70 and R2 = 0.80. Gene counts for the MSB samples were lower at 15,130 and 24,078. For the DNSB (− 20 °C storage) samples the Pearson correlations were R2 = 0.92 and R2 = 0.91. The detected gene counts for DNSB samples were 34,515 and 31,857. The Pearson correlations for all possible sample pairings in this sub-study are shown in Supplemental File 7. Scatter plots and Venn diagrams representing individual sample comparisons are in Supplemental File 8. Providing further support of our observation that storage conditions greatly influence the reproducibility of experiments, the overall Pearson correlations were lower for the MSB replicates than that for the DNSB replicates. In addition, the Venn diagrams reveal that more genespecific probes are detected in the DNSB samples compared with the MSB samples and that replicates share higher overlapping gene number
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
5
Fig. 4. Pearson correlations and averaged Venn diagrams of technical replicates from two years of three samples each.
detections in the DNSB samples compared with the MSB samples. The scatter plots confirm the Pearson correlations. Over the entire dynamic range, the replicates of gene-specific probe intensities are more reproducible for the DNSB samples compared with the MSB samples. Scatter plots (Supplemental File 8) indicate that discrepancies in the number of detected genes are concentrated near the lower detection boundaries of the array and, therefore, represent the most unreliable part of the dataset. 3.5. Sub-study 4) reducing input material The fourth sub-study explores whether reducing input material influences reproducibility. The amount of RNA extracted is proportional to the quantity of DBS sample used. Therefore, we varied parameters of DBS size and number of discs while analyzing two DNSB samples collected in the year 2000 and stored exclusively at −20 °C. The Venn diagrams in Fig. 5 illustrate the impact of varying amounts of input material on the number of detected genes. The Pearson correlations for these experiments are in Supplemental File 9. In Fig. 5, analysis of larger spot sizes produces a higher detection rate of gene-specific probes compared with smaller DBS samples. Although these data are inconclusive due to the low sample number, we infer from the Venn diagrams that “shared” expression correlates with high probe intensities (refer to the numbers in square brackets). Probe intensities are highest for gene-specific probes “shared” by all three experiments, followed by gene-specific probes shared by two experiments and then gene-specific probes unique to the individual experiments. The observation that low probe intensities in this case prompts differential gene detection between experiments is in agreement with the scatter plots (Supplemental File 8), in which differences appear concentrated near the lower detection boundaries (near the noise) of the array. Neonatal DBS samples are a limited resource intended to facilitate the analysis of DBS samples for the benefit of children and their families.
To maintain this diagnostic potential, the amount of material available for research is very limited; typically, one or two 3.2 mm spots per patient for a single project [11]. In this sub-study, varying the amounts of input material, i.e., the number and size of spots used, had minor effects (if any). No significant trend was observed in this study due to the small sample size. Any differences detected were concentrated near the lower detection boundaries of the array. Therefore, we infer that varying inputs (2 × 1.6 mm, 1 × 3.2 mm, or 2 × 3.2 mm) have a negligible impact on the analysis. Furthermore, the intensity of differentially counted genespecific probes suggests that any detection differences due to disc size would be limited because filtering often excludes probes with low intensity and little change. However, the Pearson correlations were higher in experiments using similar amounts of input material (Supplemental File 9). Therefore, maintaining a consistent setup throughout a study will improve reproducibility.
4. Conclusions DBS samples collected on filter paper, such as those collected for neonatal screenings in the Danish Neonatal Screening Biobank, are a reproducible source of RNA for expression studies. Results are significantly improved if samples are transferred to storage at −20 °C as quickly as possible following neonatal screening. We conclude that transferring DBS samples to −20 °C after neonatal screening improves the preservation of RNA. Therefore, we recommend that all similar biobanks store their samples at −20 °C. From the data presented here, we conclude that reproducible studies of whole transcriptome neonatal RNA profiles using microarrays are possible, even when samples are taken decades apart. The length of storage is secondary to the storage temperature; there were no significant differences between samples drawn in 1990 and those drawn in 2000 when all samples were stored at − 20 °C. It is also important to
Fig. 5. Venn diagrams illustrating the effect of using varying amounts of input material, i.e. spot size and number of spots, on the reproducibility. Numbers within brackets in the Venn diagrams represent the total number of detected probes in the specific sample; while numbers with no brackets are the individual number of probes in the given intersection of the Venn diagram. Numbers enclosed by square brackets constitute the mean log2 intensities of the group.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
6
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
use a single WTA kit throughout experiments because the use of multiple kits introduces significant differences in expression. Our results suggest that the amount of input material can be substantially decreased to preserve this valuable resource. Overall, this study suggests that DBS samples from the Danish Neonatal Screening Biobank can be used in research projects determining GEP via RNA microarrays. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.ymgme.2015.06.011. Authors' contributions JG was involved in the study design, analyzed and interpreted the major part of the data and drafted the manuscript. SKK extracted the RNA, conducted the array preparation and scanning, and participated in the interpretation of the results. RZN participated in the data analysis and participated in the interpretation of the results. JBP, CHSH and MBH participated in the handling of DNSB samples and the interpretation of the results. MVH and DMH initiated the study, participated in its design, the interpretation of results and drafting of the manuscript. All authors were involved in the critical revision and approved the final manuscript. Competing interests None declared. Acknowledgments The authors would like to thank and acknowledge Professor Nigel Paneth (Michigan State University) for providing access to the samples from the Michigan Screening Biobank. References [1] P. Carninci, J. Yasuda, Y. Hayashizaki, Multifaceted mammalian transcriptome, Curr. Opin. Cell Biol. 20 (3) (2008) 274–280. [2] E.D. Green, A. Chakravarti, The human genome sequence expedition: views from the “base camp”, Genome Res. 11 (5) (2001) 645–651. [3] S. Guil, M. Esteller, RNA–RNA interactions in gene regulation: the coding and noncoding players, Trends Biochem. Sci. 40 (5) (2015) 248–256. [4] Y.A. Kim, T.M. Przytycka, Bridging the gap between genotype and phenotype via network approaches, Front. Genet. 3 (2012) 227. [5] Q. Huang, Genetic study of complex diseases in the post-GWAS era, J. Genet. Genomics 42 (3) (2015) 87–98. [6] D. Welter, et al., The NHGRI GWAS Catalog, a curated resource of SNP–trait associations, Nucleic Acids Res. 42 (Database issue) (2014) D1001–D1006.
[7] Z. Dong, Y. Chen, Transcriptomics: advances and approaches, Sci. China Life Sci. 56 (10) (2013) 960–967. [8] E.L. van Dijk, et al., Ten years of next-generation sequencing technology, Trends Genet. 30 (9) (2014) 418–426. [9] J. Xuan, et al., Next-generation sequencing in the clinic: promises and challenges, Cancer Lett. 340 (2) (2013) 284–295. [10] L.C. Thygesen, et al., Introduction to Danish (nationwide) registers on health and social issues: structure, access, legislation, and archiving, Scand. J. Public Health 39 (7 Suppl.) (2011) 12–16. [11] B. Norgaard-Pedersen, D.M. Hougaard, Storage policies and use of the Danish Newborn Screening Biobank, J. Inherit. Metab. Dis. 30 (4) (2007) 530–536. [12] M.V. Hollegaard, et al., Genome-wide scans using archived neonatal dried blood spot samples, BMC Genomics 10 (2009) 297. [13] M.V. Hollegaard, et al., Robustness of genome-wide scanning using archived dried blood spot samples as a DNA source, BMC Genet. 12 (2011) 58. [14] M.V. Hollegaard, et al., DNA methylome profiling using neonatal dried blood spot samples: a proof-of-principle study, Mol. Genet. Metab. 108 (4) (2013) 225–231. [15] M.V. Hollegaard, et al., Archived neonatal dried blood spot samples can be used for accurate whole genome and exome-targeted next-generation sequencing, Mol. Genet. Metab. 110 (1–2) (2013) 65–72. [16] A.D. Borglum, et al., Genome-wide study of association and interaction with maternal cytomegalovirus infection suggests new schizophrenia loci, Mol. Psychiatry 19 (3) (2014) 325–333. [17] M. Horikoshi, et al., New loci associated with birth weight identify genetic links between intrauterine growth and adult height and metabolism, Nat. Genet. 45 (1) (2013) 76–82. [18] S. Steinberg, et al., Common variant at 16p11.2 conferring risk of psychosis, Mol. Psychiatry 19 (1) (2014) 108–114. [19] K. Bonnelykke, et al., A genome-wide association study identifies CDHR3 as a susceptibility locus for early childhood asthma with severe exacerbations, Nat. Genet. 46 (1) (2014) 51–55. [20] F. Geller, et al., Genome-wide association analyses identify variants in developmental genes associated with hypospadias, Nat. Genet. 46 (9) (2014) 957–963. [21] B. Feenstra, et al., Common variants associated with general and MMR vaccinerelated febrile seizures, Nat. Genet. 46 (12) (2014) 1274–1282. [22] Y.H. Zhang, E.R. McCabe, RNA analysis from newborn screening dried blood specimens, Hum. Genet. 89 (3) (1992) 311–314. [23] S.K. Khoo, et al., Acquiring genome-wide gene expression profiles in Guthrie card blood spots using microarrays, Pathol. Int. 61 (1) (2011) 1–6. [24] F. Gauffin, et al., Quantitation of RNA decay in dried blood spots during 20 years of storage, Clin. Chem. Lab. Med. 47 (12) (2009) 1467–1469. [25] J.H. Resau, et al., Evaluation of sex-specific gene expression in archived dried blood spots (DBS), Int. J. Mol. Sci. 13 (8) (2012) 9599–9608. [26] C. Wei, et al., Comparison of frozen and unfrozen blood spots for gene expression studies, J. Pediatr. 164 (1) (2014) 189–191.e1. [27] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate — a practical and powerful approach to multiple testing, J. R. Stat. Soc. Ser. B Methodol. 57 (1) (1995) 289–300. [28] Y. Benjamini, D. Yekutieli, The Control of the False Discovery Rate in Multiple Testing Under Dependency, 2001. 1165–1188. [29] J. Hedegaard, et al., Next-generation sequencing of RNA and DNA isolated from paired fresh-frozen and formalin-fixed paraffin-embedded samples of human cancer and normal tissue, PLoS One 9 (5) (2014) e98187. [30] A.M. Morse, et al., Comparison between NuGEN's WT-Ovation Pico and one-direct amplification systems, J. Biomol. Tech. 21 (3) (2010) 141–147.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
Contents lists available at ScienceDirect
Molecular Genetics and Metabolism journal homepage: www.elsevier.com/locate/ymgme
Regular Article
Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank Jonas Grauholm a, Sok Kean Khoo b, Radoslav Z. Nickolov c, Jesper B. Poulsen a, Marie Bækvad-Hansen a, Christine S. Hansen a, David M. Hougaard d, Mads V. Hollegaard a,⁎ a
Section of Neonatal Genetics, Danish Centre for Neonatal Screening, Department of Congenital Diseases, Statens Serum Institut, Artillerivej 5, DK-2300 Copenhagen S, Denmark Department of Cell and Molecular Biology, Grand Valley State University, Grand Rapids, MI 49503, USA Department of Mathematics & Computer Science, Fayetteville State University, Fayetteville, NC 28301, USA d Danish Centre for Neonatal Screening, Department of Congenital Diseases, Statens Serum Institut, Artillerivej 5, DK-2300 Copenhagen S, Denmark b c
a r t i c l e
i n f o
Article history: Received 5 May 2015 Received in revised form 26 June 2015 Accepted 26 June 2015 Available online xxxx Keywords: Biobank Dried blood spots Neonatal screening RNA storage Gene expression profile Whole transcriptome amplification
a b s t r a c t A large part of the human genome is transcribed into various forms of RNA, and the global gene expression profile (GEP) has been studied for several years using technology such as RNA-microarrays. In this study, we evaluate whether neonatal dried blood spot (DBS) samples stored in the Danish Neonatal Screening Biobank (DNSB) can be used for GEP. This paper is divided into sub-studies examining the effects of: 1) different whole transcriptome amplification kits (WTA); 2) years of storage and storage in room temperature (RT) versus freezers (−20 °C) on DNSB DBS samples; 3) effects of RT storage vs freezer storage on DBS samples from the USA and DNSB, and 4) using smaller disc sizes, thereby decreasing DBS use. We present evidence that reliable and reproducible GEPs can be obtained using neonatal DBS samples. The main source of variation is the storage condition. When samples are stored at −20 °C, the dynamic range is increased, and Pearson correlations are higher. Differential analysis reveals no statistically significant differences between samples collected a decade apart and stored at −20 °C. However, samples stored at RT show differential expression for a third of the gene-specific probes. Our data also suggests that using alternate WTA kits significantly changes the GEP. Finally, the amount of input material, i.e., the size and number of DBS discs used, can be reduced to preserve this valuable and limited material. We conclude that DNSB DBS samples provide a reproducible resource for GEP. Results are improved if the cards are stored at −20 °C. Furthermore, it is important to use a single type of kit for analysis because using alternate kits introduces differential expression. © 2015 Elsevier Inc. All rights reserved.
1. Introduction RNA is a diverse molecule essential to cellular function. One form carries the genetic code from the nucleus while others control diverse regulatory functions, such as splicing and translation. The ENCODE project showed that up to 93% of the human genome can be transcribed into one of several known forms of RNA [1–3]. Since the introduction of genome wide association studies, approximately 2000 studies have linked numerous diseases to more than 10,000 SNPs. The challenge now is to understand their functional consequences. Especially in complex diseases, it is hard to infer organism level ⁎ Corresponding author at: Section of Neonatal Genetics, Danish Center for Neonatal Screening, Department of Congenital Diseases, Statens Serum Institut, Artillerivej 5, DK-2300 Copenhagen S, Denmark. E-mail addresses: [email protected] (J. Grauholm), [email protected] (S.K. Khoo), [email protected] (R.Z. Nickolov), [email protected] (J.B. Poulsen), [email protected] (M. Bækvad-Hansen), [email protected] (C.S. Hansen), [email protected] (D.M. Hougaard), [email protected] (M.V. Hollegaard).
phenotypes (i.e., disease state) from observed genotypes. Changes in molecular phenotype, as observed through up- or down-regulation of RNA, may help bridge the gap between genotypes and phenotypes [4–6]. Numerous technologies and methodologies are used to study RNA, and studies of gene expression profiles (GEP) are expected to have an increasing impact in research and as a clinical tool [7–9]. Denmark has a long tradition of collecting and storing data in nationwide registries. The registries, which can be linked through a unique person identifier number (CPR-number), contain information from birth to death, including notable social and health-related events [10]. The accessibility of well-defined samples is generally a limiting factor in genetics. In addition to the data registries, Denmark also has a long history of storing biological material, including the excess neonatal dried blood spot (DBS) samples routinely collected as a part of the Danish Neonatal Screening program. The DBS samples, stored in the Danish Neonatal Screening Biobank (DNSB), are primarily used to improve neonatal screening and diagnostics [11]. It is also possible to access the DBS samples for research on improving the Danish population's
http://dx.doi.org/10.1016/j.ymgme.2015.06.011 1096-7192/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
2
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
health in general. By January 2015, the DNSB contained more than 2.2 million DBS samples dating back to 1982. Combining the biological material from the DNSB DBS samples with data from the registries unlocks a powerful tool for conducting unbiased population-based health research which are based on well-documented diagnoses and thorough anamnesis. We have previously shown that DNSB DBS samples can be used for array genotyping [12,13], DNA methylation [14], and whole genome and exome sequencing [15]. Using these approaches, the DNSB has been successfully used to investigate several diseases [16–21]. Other groups have published RNA work using DBS cards as a source [22–24]. A series of work by Dr. Khoo and colleagues has convincingly shown that DBS samples from the Michigan Screening Biobank (MSB) can be used for GEP [23] and to detect the differential expression of gender specific genes [25]. Dr. Khoo and her colleagues have also documented differences between DBS stored at room temperature versus conventionally preserved blood after more than 10 years in storage [26]. This study evaluates whether DBS samples from the DNSB can be used for reliable GEPs. Moreover, we evaluate the effect of storage conditions, including time and temperature, the effect of using different whole-transcriptome amplification (WTA) kits and whether the GEP remains accurate and stable when using smaller disc sizes, which provide less whole blood for the RNA extraction.
2. Materials and methods
Denmark. The Institutional Review Board at the Michigan State University approved the use of neonatal DBS samples from the MSB. This study consists of four sub-units (see Fig. 1) evaluating the effect on GEPs when: 1) Using two different WTA kits to compare samples collected in the year 2000 and stored at −20 °C. 2) Testing DBS samples stored in different conditions, using two sets of DNSB samples collected in 1982 and 1990. While the 1990 samples were stored at −20 °C shortly after the neonatal screening analysis, the 1982 samples were stored at room temperature (RT) for up to a month before being transferred to −20 °C. Because samples collected in 2000 were stored in the same conditions as the 1990 samples, we also compare the results of this analysis to data generated in Sub-study 1 to determine whether the duration of storage causes variance. 3) Comparing DBS samples stored at RT to those stored at −20 °C. To provide further support that storage conditions, rather than years of storage, are the primary source of variance, we compare technical replicates of neonatal DBS samples collected in 1984 and 2011 and stored at RT in the MSB. As a reference, we use two DNSB neonatal DBS samples collected in 2000 and stored exclusively at −20 °C. 4) Decreasing the amount of whole-blood (WB) used for RNA extraction and, thereby, the WTA RNA input. We analyzed DNSB DBS samples stored at − 20 °C since the year 2000. Three conditions were tested: 2 × 3.2 mm discs (~ 6 μL WB), 1 × 3.2 mm discs (~3 μL WB), and 2 × 1.6 mm discs (~1.5 μL WB).
2.1. Sample material and study design 2.2. RNA isolation, WTA, and microarray technology In this study, we used DBS samples stripped of all information except the year of storage. As such, this study does not constitute a healthrelated research project as defined by the Danish “Act on Research Ethics Review of Health Research Projects”, but is considered a developmental project for the newborn screening program. Development projects for neonatal screening do not require prior approval from the Committees on Biomedical Research Ethics for the Capital Region of
Total RNA was extracted, purified, and concentrated using an Illustra RNAspin mini kit (GE Healthcare, Little Chalfont, Buckinghamshire, United Kingdom) and treated with DNase to eliminate DNA contamination. We used the WT-Ovation Pico RNA amplification system (NuGEN) (NuGEN, San Carlos, California, USA) or the Sigma-Aldrich's complete whole transcriptome amplification kit (Sigma) (Sigma, St. Louis,
Fig. 1. Schematic overview of storage conditions and sample processing of the four sub-units in this study. The MSB neonatal DBS samples were stored at RT prior to the experiment (2 and 30 years respectively). The DNSB neonatal DBS samples are stored predominantly at −20 °C since 1982, 1990 and 2000 (12, 20 and 30 years, respectively). The DBS samples were excised in sizes of 2 × 3.2 mm, 1 × 3.2 mm or 2 × 1.6 mm, their RNA extracted and WTA using NuGEN or Sigma kits, before subjection of the samples to microarray analysis with the Agilent 8 × 60K RNA array. Sub-study 1 compares the effects of WTA kits. Sub-study 2 compares the storage conditions used in the DNSB. Sub-study 3 compare DNSB samples stored at frost with MSB samples stored at RT. Sub-study 4 explores if it is possible to reduce the amount of input material.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
3
Missouri, USA) to generate single-stranded cDNA that was subsequently labeled with the fluorescent dye Alexa Fluor 3 before hybridization onto the Agilent whole human genome gene expression 8 × 60K microarray (Agilent 60K) (Agilent, Santa Clara, California, USA). Each array contains 62,976 oligonucleotide probes covering 27,958 Entrez Gene RNAs and 7419 long intergenic non-coding RNAs. The arrays were hybridized for 17 h at 65 °C and a rotation speed of 10 rpm, washed for 2 min each with wash buffers 1 and 2 (Agilent, Santa Clara, California, USA), and scanned using an Agilent G3 high-resolution scanner (Agilent, Santa Clara, California, USA). Probe features were extracted from the microarray data using Agilent's Feature Extraction software v.10.7.3.1 (Agilent), and used for GEP analysis via Partek Genomics Suite (PGS) version 6.6 (Partek Incorporated, St. Louis, Missouri, USA). 2.3. Data analysis The gProcessedSignal and gIsInNegCtrlRange values were loaded into PGS. The gProcessedSignal measurement is an analogue value based on a given probe's intensity and is assumed equivalent to the expression level of the transcript corresponding to the probe. The gIsInNegCtrlRange measurement is a digital value that is set as “1” when the intensity value of the gProcessedSignal is within three standard deviations of the array background. Probes that were not flagged (i.e., those with gIsInNegCtrlRange values of 0) are considered detected. The Initial QC was based on the gProcessedSignal value for the unfiltered set of ~60K probes. Evaluation through principal component analysis (PCA) was based on a covariance model and hierarchical clustering analysis (HCA) based on a Euclidean model of complete linkage. To generate Venn diagrams, we quantified the shared and unique entries in the lists of detected probes. We then log2-transformed and quantile-normalized the gProcessedSignal, which was subsequently used to generate intensity–intensity scatterplots of sample correlations and to perform an analysis of variance (ANOVA). The ANOVAs were run as separate one-way analyses, where the factors were the age of the card or the WTA kit. Significant hits were those with at least a 2-fold expression change and a p-value of less than 0.01 following a Benjamini and Hochberg correction for false discovery rate [27,28]. 3. Results and discussion 3.1. RNA microarray quality control The first step in our RNA array quality control analysis was to consult the un-normalized log2 transformed intensities to identify potential biases and samples to exclude. The raw probe intensity QC showed an increased dynamic range for gene-specific probes with higher average intensities compared with the negative control probes (Supplemental File 1); demonstrating that gene-specific probe signals are above background levels in all samples. Not only does this confirm the presence of reliable material in all samples used in the microarray analysis, but it also validates the technical aspects of the experimental setup. Following quantile normalization of the intensity values, we used PGS to perform a PCA using a covariance model. Using the normalized data in PGS, we created an HCA using a Euclidean model of complete linkage. For both methods, we used the complete set of gene-specific probes to avoid selection bias. The PCA and HCA (Fig. 2) show distinctive clustering of the groups. Based on the PCA (Fig. 2A), we were able to define the two main contributors of variance. The storage condition (i.e., − 20 °C versus RT) of DBS samples accounts for 51.9% of the variation (PC1) and the use of different WTA kits accounts for 13% of the variation (PC2). Combined, the two principal components separate samples into groups of four quadrates with non-overlapping separation. The quadrates, separated by the midpoint of the axes, only contain samples from one of the four possible groupings based on the analysis with two WTA kits combined and two storage conditions.
Fig. 2. A) PCA generated using a covariance model. B) HCA constructed using a Euclidian model for complete linkage. Both (A & B) are generated using all gene-specific probes as to avoid selection bias. Each dot in the PCA plot and each line in the HC represent one sample, individual sample names have been removed, instead we have color coded by storage and preparation kit: Sigma-(−20) (green), Sigma-RT (purple), NuGEN-(−20) (red) and NuGEN-RT (blue). Principal component 1 (PC1) stratifies the samples on storage conditions (−20 °C vs RT), whereas PC2 splits by the kit used on the WTA (Sigma vs NuGEN). The HCA shows a similar dependency by WTA kit and storage conditions as in the PCA plot, though the samples stored at room temperature (RT) group as a single branch in the HCA.
Similar to the PCA, the HCA in Fig. 2B shows sample dependency by WTA kit and storage condition. However, the RT storage samples cluster (Fig. 2B). Importantly, when repeating the analysis using only Agilent control probes, neither the PCA nor the HCA produces the pattern shown in Fig. 2 (see Supplemental File 2). This pattern difference suggests that the effect seen in Fig. 2 is the result of variance from storage conditions and the WTA kit used, rather than artifacts introduced by the experimental setup. Our observations are consistent with those previously reported. Storage conditions have been established as a major factor in RNA preservation e.g., in the extreme case, comparing formalin-fixed, paraffinembedded samples to their fresh frozen counterparts [29]. Similarly, variances from the use of different WTA kits are well described [30]. In contrast, the duration of DBS sample storage did not play an important role in the observed variances. In both the PCA and HCA, the samples drawn in 1990 cluster with the samples drawn in 2000. Thus, we conclude that storage condition, rather than storage duration, is the primary source of high variance in DBS sample analysis. As a result, we are confident that our measurements reflect the biology of samples rather than an introduced bias; we therefore analyzed the full data set. With respect to storage condition, differences appear between groups, even at the raw data level. In Supplemental File 1, the signalto-noise ratio increases in samples from − 20 °C versus RT storage. Extracting the median signal for all RT or − 20 °C negative controls, we found a probe intensity just above 2 (in both), while the gene-specific probes had average median signals of 3.7 in the RT and 4.8 in the −20 °C, confirming that the signal to noise ratio is increased by storage at −20 °C.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
4
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
3.2. Sub-study 1) whole transcriptome amplification kits Eight samples were analyzed as part of the first sub-study. These were all DNSB samples from the year 2000 that were stored exclusively at −20 °C. Transcription profiles were generated for each sample based on the use of two WTA kits. Consequently, identical samples were processed twice, once using the NuGEN kit and once using the Sigma kit. For each sample, we quantified gene-specific probes unique to one of the kits and those detected in both. The counts were averaged between samples to create a single Venn diagram for all comparisons. We also calculated the Pearson correlations, listed with the Venn diagram in Fig. 3. Fig. 3 shows the cross-kit replications comparing identical samples. The Pearson correlations for all combinations of WTA kit and samples are in Supplemental File 3. Overall, the Pearson correlation is higher between different samples tested with the same WTA kit (Pearson average correlation ± standard deviation for NuGEN: R2 = 0.95 ± 0.00 and Sigma: R2 = 0.97 ± 0.00) than between identical samples amplified using different WTA kits (R2 = 0.91 ± 0.01). The raw counts of unique and commonly detected genes used to calculate the mean values in the Venn diagram (Fig. 3) are shown in Supplemental File 4. We consistently detected a higher number of genes using the Sigma kit versus the NuGEN kit. Next, we did an ANOVA to compare identical sample sets prepared using different WTA kits. After having applied a Benjamini and Hochberg correction [27,28] for false discovery rate, we identified 8462 genes with at least a 2-fold change in expression level and a significance of P b 0.01. These data provide further support that, for experimental groups to remain comparable with respect to specific probe detection and expression levels, it is essential to prepare samples using the same WTA kit.
3.3. Sub-study 2) effects of storage time and temperature The second sub-study further investigates the reproducibility of technical replicates and explores the effect of storage at −20 °C versus RT. This analysis used two sets of DNSB samples: the first collected in 1982 and the second in 1990. Samples collected in 1982 had a latency of up to a month before being transferred to − 20 °C, whereas those collected in 1990 were transferred to −20 °C shortly after screening. We quantified the gene-specific probes detected as unique to one or common to both of the replicates. These were averaged to create a single Venn diagram, which is shown in Fig. 4 along with the Pearson correlations. The Pearson correlations comparing all samples in the sub-study are available in Supplemental File 5. The raw counts of the unique and commonly detected genes used to calculate the means shown in the Venn diagrams in Fig. 4 are in Supplemental File 6.
Fig. 3. Pearson correlations and an averaged Venn diagram for the samples used to compare the WTA kits.
Of the two groups with higher Pearson correlation of technical replicates, there is a significant difference between the 1990 samples (R2 = 0.97 ± 0.01) and the 1982 samples (R2 = 0.91 ± 0.02). Furthermore, Supplemental File 5 shows that the Pearson correlation between any two random samples in the 1990 cohort (R2 = 0.97 ± 0.01) is better than the correlation between technical replicates from the 1982 sample mentioned above. Although we did not study technical replicates for samples from the year 2000, we again observed improved Pearson correlations between unrelated individual samples stored at − 20 °C. In Supplemental File 3, R2 = 0.97 ± 0.00 when comparing any two unrelated samples prepared using the Sigma kit. These data indicate that performance is more consistent for samples directly transferred to storage − 20 °C, i.e., the 1990 and 2000 samples, despite the 10year difference in storage duration. Thus, samples whose collection is separated by decades can be included in the same study as long as storage and experimental conditions are stringently controlled. There is however, a marked reduction in reproducibility when analyzing the 1982 samples, which had been stored at room temperature for up to a month before being transferred to storage at −20 °C for three decades. Pearson correlations for the 1982 sample group are in accordance with previous reports on samples stored at room temperature for longer periods. For similar samples tested on the Agilent 44 k array, Khoo et al. found a correlation for experimental replicates of R2 = 0.92 [23]. Overall, the values reported here, R2 = 0.97 ± 0.01 and R2 = 0.91 ± 0.02 for different sample groups, are above those of experimental replicates reported as satisfactory for FFPE samples (R2 = 0.9) [29]. To further investigate the effect of storage time, we performed an ANOVA analysis on DNSB samples from the years 1982, 1990 and 2000. We did not identify any genes with differential expression while comparing DNSB samples from 1990 with those from 2000. In contrast, we detected 22,356 differentially expressed gene-specific probes of the 58,782 probes tested when comparing samples from 1982 to 1990. A comparison of samples from 1982 and those from 2000 detected 27,660 differentially expressed gene-specific probes. For all comparisons, hits were only retained if the p-value was less than 0.01 following Benjamini and Hochberg correction and if the hit had at least a two-fold change in signal intensity between samples (i.e., a log2 change of ±1). These experiments demonstrate that the RNA in DBS samples stored at −20 °C serves is better preserved, leading to increased success and accuracy in potential RNA array studies. 3.4. Sub-study 3) storage at RT versus −20 °C To further probe the effect of room temperature storage on the reliability of DBS samples used in RNA profiling, we carried out experiments using two MSB samples that had been stored exclusively at room temperatures. To highlight the differences between RT and −20 °C storage, we compared the MSB samples to two DNSB samples. Samples were prepared using the Sigma and the NuGEN kit and, while we have already shown reduced decreased performance in cross kit comparisons from the PCA (see Fig. 2), we expect the storage conditions to be a bigger cause of variation. The Pearson correlations for the MSB samples (RT storage) were R2 = 0.70 and R2 = 0.80. Gene counts for the MSB samples were lower at 15,130 and 24,078. For the DNSB (− 20 °C storage) samples the Pearson correlations were R2 = 0.92 and R2 = 0.91. The detected gene counts for DNSB samples were 34,515 and 31,857. The Pearson correlations for all possible sample pairings in this sub-study are shown in Supplemental File 7. Scatter plots and Venn diagrams representing individual sample comparisons are in Supplemental File 8. Providing further support of our observation that storage conditions greatly influence the reproducibility of experiments, the overall Pearson correlations were lower for the MSB replicates than that for the DNSB replicates. In addition, the Venn diagrams reveal that more genespecific probes are detected in the DNSB samples compared with the MSB samples and that replicates share higher overlapping gene number
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
5
Fig. 4. Pearson correlations and averaged Venn diagrams of technical replicates from two years of three samples each.
detections in the DNSB samples compared with the MSB samples. The scatter plots confirm the Pearson correlations. Over the entire dynamic range, the replicates of gene-specific probe intensities are more reproducible for the DNSB samples compared with the MSB samples. Scatter plots (Supplemental File 8) indicate that discrepancies in the number of detected genes are concentrated near the lower detection boundaries of the array and, therefore, represent the most unreliable part of the dataset. 3.5. Sub-study 4) reducing input material The fourth sub-study explores whether reducing input material influences reproducibility. The amount of RNA extracted is proportional to the quantity of DBS sample used. Therefore, we varied parameters of DBS size and number of discs while analyzing two DNSB samples collected in the year 2000 and stored exclusively at −20 °C. The Venn diagrams in Fig. 5 illustrate the impact of varying amounts of input material on the number of detected genes. The Pearson correlations for these experiments are in Supplemental File 9. In Fig. 5, analysis of larger spot sizes produces a higher detection rate of gene-specific probes compared with smaller DBS samples. Although these data are inconclusive due to the low sample number, we infer from the Venn diagrams that “shared” expression correlates with high probe intensities (refer to the numbers in square brackets). Probe intensities are highest for gene-specific probes “shared” by all three experiments, followed by gene-specific probes shared by two experiments and then gene-specific probes unique to the individual experiments. The observation that low probe intensities in this case prompts differential gene detection between experiments is in agreement with the scatter plots (Supplemental File 8), in which differences appear concentrated near the lower detection boundaries (near the noise) of the array. Neonatal DBS samples are a limited resource intended to facilitate the analysis of DBS samples for the benefit of children and their families.
To maintain this diagnostic potential, the amount of material available for research is very limited; typically, one or two 3.2 mm spots per patient for a single project [11]. In this sub-study, varying the amounts of input material, i.e., the number and size of spots used, had minor effects (if any). No significant trend was observed in this study due to the small sample size. Any differences detected were concentrated near the lower detection boundaries of the array. Therefore, we infer that varying inputs (2 × 1.6 mm, 1 × 3.2 mm, or 2 × 3.2 mm) have a negligible impact on the analysis. Furthermore, the intensity of differentially counted genespecific probes suggests that any detection differences due to disc size would be limited because filtering often excludes probes with low intensity and little change. However, the Pearson correlations were higher in experiments using similar amounts of input material (Supplemental File 9). Therefore, maintaining a consistent setup throughout a study will improve reproducibility.
4. Conclusions DBS samples collected on filter paper, such as those collected for neonatal screenings in the Danish Neonatal Screening Biobank, are a reproducible source of RNA for expression studies. Results are significantly improved if samples are transferred to storage at −20 °C as quickly as possible following neonatal screening. We conclude that transferring DBS samples to −20 °C after neonatal screening improves the preservation of RNA. Therefore, we recommend that all similar biobanks store their samples at −20 °C. From the data presented here, we conclude that reproducible studies of whole transcriptome neonatal RNA profiles using microarrays are possible, even when samples are taken decades apart. The length of storage is secondary to the storage temperature; there were no significant differences between samples drawn in 1990 and those drawn in 2000 when all samples were stored at − 20 °C. It is also important to
Fig. 5. Venn diagrams illustrating the effect of using varying amounts of input material, i.e. spot size and number of spots, on the reproducibility. Numbers within brackets in the Venn diagrams represent the total number of detected probes in the specific sample; while numbers with no brackets are the individual number of probes in the given intersection of the Venn diagram. Numbers enclosed by square brackets constitute the mean log2 intensities of the group.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011
6
J. Grauholm et al. / Molecular Genetics and Metabolism xxx (2015) xxx–xxx
use a single WTA kit throughout experiments because the use of multiple kits introduces significant differences in expression. Our results suggest that the amount of input material can be substantially decreased to preserve this valuable resource. Overall, this study suggests that DBS samples from the Danish Neonatal Screening Biobank can be used in research projects determining GEP via RNA microarrays. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.ymgme.2015.06.011. Authors' contributions JG was involved in the study design, analyzed and interpreted the major part of the data and drafted the manuscript. SKK extracted the RNA, conducted the array preparation and scanning, and participated in the interpretation of the results. RZN participated in the data analysis and participated in the interpretation of the results. JBP, CHSH and MBH participated in the handling of DNSB samples and the interpretation of the results. MVH and DMH initiated the study, participated in its design, the interpretation of results and drafting of the manuscript. All authors were involved in the critical revision and approved the final manuscript. Competing interests None declared. Acknowledgments The authors would like to thank and acknowledge Professor Nigel Paneth (Michigan State University) for providing access to the samples from the Michigan Screening Biobank. References [1] P. Carninci, J. Yasuda, Y. Hayashizaki, Multifaceted mammalian transcriptome, Curr. Opin. Cell Biol. 20 (3) (2008) 274–280. [2] E.D. Green, A. Chakravarti, The human genome sequence expedition: views from the “base camp”, Genome Res. 11 (5) (2001) 645–651. [3] S. Guil, M. Esteller, RNA–RNA interactions in gene regulation: the coding and noncoding players, Trends Biochem. Sci. 40 (5) (2015) 248–256. [4] Y.A. Kim, T.M. Przytycka, Bridging the gap between genotype and phenotype via network approaches, Front. Genet. 3 (2012) 227. [5] Q. Huang, Genetic study of complex diseases in the post-GWAS era, J. Genet. Genomics 42 (3) (2015) 87–98. [6] D. Welter, et al., The NHGRI GWAS Catalog, a curated resource of SNP–trait associations, Nucleic Acids Res. 42 (Database issue) (2014) D1001–D1006.
[7] Z. Dong, Y. Chen, Transcriptomics: advances and approaches, Sci. China Life Sci. 56 (10) (2013) 960–967. [8] E.L. van Dijk, et al., Ten years of next-generation sequencing technology, Trends Genet. 30 (9) (2014) 418–426. [9] J. Xuan, et al., Next-generation sequencing in the clinic: promises and challenges, Cancer Lett. 340 (2) (2013) 284–295. [10] L.C. Thygesen, et al., Introduction to Danish (nationwide) registers on health and social issues: structure, access, legislation, and archiving, Scand. J. Public Health 39 (7 Suppl.) (2011) 12–16. [11] B. Norgaard-Pedersen, D.M. Hougaard, Storage policies and use of the Danish Newborn Screening Biobank, J. Inherit. Metab. Dis. 30 (4) (2007) 530–536. [12] M.V. Hollegaard, et al., Genome-wide scans using archived neonatal dried blood spot samples, BMC Genomics 10 (2009) 297. [13] M.V. Hollegaard, et al., Robustness of genome-wide scanning using archived dried blood spot samples as a DNA source, BMC Genet. 12 (2011) 58. [14] M.V. Hollegaard, et al., DNA methylome profiling using neonatal dried blood spot samples: a proof-of-principle study, Mol. Genet. Metab. 108 (4) (2013) 225–231. [15] M.V. Hollegaard, et al., Archived neonatal dried blood spot samples can be used for accurate whole genome and exome-targeted next-generation sequencing, Mol. Genet. Metab. 110 (1–2) (2013) 65–72. [16] A.D. Borglum, et al., Genome-wide study of association and interaction with maternal cytomegalovirus infection suggests new schizophrenia loci, Mol. Psychiatry 19 (3) (2014) 325–333. [17] M. Horikoshi, et al., New loci associated with birth weight identify genetic links between intrauterine growth and adult height and metabolism, Nat. Genet. 45 (1) (2013) 76–82. [18] S. Steinberg, et al., Common variant at 16p11.2 conferring risk of psychosis, Mol. Psychiatry 19 (1) (2014) 108–114. [19] K. Bonnelykke, et al., A genome-wide association study identifies CDHR3 as a susceptibility locus for early childhood asthma with severe exacerbations, Nat. Genet. 46 (1) (2014) 51–55. [20] F. Geller, et al., Genome-wide association analyses identify variants in developmental genes associated with hypospadias, Nat. Genet. 46 (9) (2014) 957–963. [21] B. Feenstra, et al., Common variants associated with general and MMR vaccinerelated febrile seizures, Nat. Genet. 46 (12) (2014) 1274–1282. [22] Y.H. Zhang, E.R. McCabe, RNA analysis from newborn screening dried blood specimens, Hum. Genet. 89 (3) (1992) 311–314. [23] S.K. Khoo, et al., Acquiring genome-wide gene expression profiles in Guthrie card blood spots using microarrays, Pathol. Int. 61 (1) (2011) 1–6. [24] F. Gauffin, et al., Quantitation of RNA decay in dried blood spots during 20 years of storage, Clin. Chem. Lab. Med. 47 (12) (2009) 1467–1469. [25] J.H. Resau, et al., Evaluation of sex-specific gene expression in archived dried blood spots (DBS), Int. J. Mol. Sci. 13 (8) (2012) 9599–9608. [26] C. Wei, et al., Comparison of frozen and unfrozen blood spots for gene expression studies, J. Pediatr. 164 (1) (2014) 189–191.e1. [27] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate — a practical and powerful approach to multiple testing, J. R. Stat. Soc. Ser. B Methodol. 57 (1) (1995) 289–300. [28] Y. Benjamini, D. Yekutieli, The Control of the False Discovery Rate in Multiple Testing Under Dependency, 2001. 1165–1188. [29] J. Hedegaard, et al., Next-generation sequencing of RNA and DNA isolated from paired fresh-frozen and formalin-fixed paraffin-embedded samples of human cancer and normal tissue, PLoS One 9 (5) (2014) e98187. [30] A.M. Morse, et al., Comparison between NuGEN's WT-Ovation Pico and one-direct amplification systems, J. Biomol. Tech. 21 (3) (2010) 141–147.
Please cite this article as: J. Grauholm, et al., Gene expression profiling of archived dried blood spot samples from the Danish Neonatal Screening Biobank, Mol. Genet. Metab. (2015), http://dx.doi.org/10.1016/j.ymgme.2015.06.011