Importance of experimental information (metadata) for archived sequence data: case of specific gene bias due to lag time between sample harvest and RNA protection in RNA sequencing

Introduction

High-throughput sequencing (HTS) has been applied in several fields of biological research since its appearance about a decade ago and has made it possible to collect data with substantially higher-throughput and lower-cost (Kircher & Kelso, 2010; Kodama, Shumway & Leinonen, 2012). Large volumes of HTS data have been submitted to the Sequence Read Archive (SRA, Leinonen, Sugawara & Shumway, 2010), and repository users can access the archived data for reuse. Nakazato, Ohta & Bono (2013) constructed a publication list including PubMed IDs (PMIDs) referring to SRA entries. This list is useful for retrieving archived data based on the experimental details (metadata); however, metadata described by the data depositor lack some important information about protocol steps (Alnasir & Shanahan, 2015). Ohta, Nakazato & Bono (2017) calculated the quality of all archived sequence data in the SRA to allow users to control not only metadata but also the data quality in their searches.

It may not be possible to clarify the validity of the sequence data by examining the quality control data generated by a quality control tool (e.g., FastQC) without metadata containing experiment-specific information. For example, this author has obtained two RNA samples: one isolated from cultured cells in 30 min after sampling, and the other isolated from the same cells in 60 min after sampling; both of these RNA samples prepared in the same laboratory showed high RNA integrity number (RIN) values indicating non-degraded RNA. However, it remains questionable whether data obtained from RNA samples with different time lags from sample harvest to RNA protection can be handled equally. In the case of RNA sequencing (RNA-Seq), it is necessary to consider that changes in gene expression can occur over a short time period (several seconds to a few minutes) in many organisms. For instance, in the filamentous fungus Aspergillus fumigatus, 23% and 35% of genes were differentially expressed after 30 and 120 min of hypoxia exposure, respectively (Losada et al., 2014) compared with before exposure. In the case of cultured cells, dynamic changes of gene expression profiles occurred during the first 2 h (30, 60, 90, and 120 min) of serum stimulation in hTERT immortalized foreskin-derived human fibroblasts (Kirkconnell et al., 2016). Therefore, to isolate RNA that accurately reflects the transcriptome at the point of harvest, raw biological samples should be frozen in liquid nitrogen or treated with RNA stabilization reagent or RNA lysis buffer as soon as possible.

Some cultured mammalian cells need to be suspended in 1 × PBS to remove fetal bovine serum and/or media before treatment with RNA lysis buffer. For example, the strategy for RNA sample processing from suspended Chinese hamster ovary (CHO) cells at the author’s laboratory at Nihon BioData Corporation is as follows: (1) The medium containing the suspended CHO cells is transferred from the cell culture flask to a 1.5 mL centrifuge tube. (2) The centrifuge tube is centrifuged for 5 min at 300×g and all supernatant is removed. (3) 200 µL of 1 × PBS is added to the cell pellet and pipetted up and down five times to thoroughly suspend the pellet. (4) The centrifuge tube is centrifuged for 5 min at 300×g and all supernatant is removed. (5) The centrifuge tube is flicked to loosen the cell pellet from the bottom. (6) 350 µL of buffer RLT (RNA lysis buffer containing guanidine thiocyanate; QIAGEN, Hilden, Germany) is added to the pellet and pipetted up and down five times to mix thoroughly. (7) The entire cell suspension is immediately loaded onto a QIAshredder spin column (QIAGEN) and homogenized by centrifugation for 2 min in a microcentrifuge at maximum speed. When cells are treated by using the above strategy, the time lag until the addition of RNA lysis buffer, which contains guanidine hydrochloride, is critical. Because RNA in cells suspended in RNA lysis buffer immediately inactivates RNases to ensure isolation of intact RNA, the time lag until the addition of this reagent should be as short as possible. If the number of samples is less than 10, RNA lysis buffer can be added to the cell pellet within about 15 min. As the number of samples increases, however, the time it takes for centrifuging and pipette operation for samples becomes longer. As a result, the cells remain in media or 1 × PBS at the conditions (temperature and CO₂ concentration) outside the incubator for a long time.

Here, to evaluate whether RNA-Seq data obtained from RNA samples with a different lag time until RNA protection can be analyzed equally, we harvested CHO-S cells after 3, 5, 6, and 7 days of cultivation, added RNA lysis buffer at various lag times (Fig. 1; 15, 30, 45, and 60 min), purified the RNA, and then performed RNA-Seq. We found that the global trends of gene expression were similar, but the expression of some specific genes was significantly different between samples that were processed with different lag times until RNA protection on the same cultivation day. Most of the genes whose profiles changed with lag time until RNA protection were related to apoptosis. Based on our results, we propose to describe detailed protocol procedures including working time for reusing sequence data accessed from the data repository.

Materials & Methods

Cells

Bach CHO-S (Thermo Fisher Scientific, Waltham, MA, USA) cell culture was grown in 20 ml of Gibco CD-CHO medium (Thermo Fisher Scientific, Waltham, MA, USA) containing 400 µL of 200 mM L-Alanyl-L-glutamine (final concentration: 4 mM; Sigma–Aldrich, St. Louis, MI, USA) and 40 µL of Gibco Anti-Clumping Agent (final concentration, 0.2% (v/v); Thermo Fisher Scientific, Waltham, MA, USA) in a 125 mL flask. The flask was shaken (120 rpm) at 37 °C in a humidified 5% CO₂ atmosphere. Determination of total and viable cell numbers was performed with Vi-Cell XR Cell Viability Analyzer (Beckman Coulter, Fullerton, CA, USA) and CellProfiler software (Carpenter et al., 2006). The images obtained from Vi-Cell XR were analyzed and total and viable cell numbers were counted using the CellProfiler pipeline. The growth curve of CHO-S cells over the period of study is shown in Fig. 2.

Results

To determine whether RNA-Seq data obtained from RNA samples that differ in the time lag from harvesting cells to protecting RNA (hereafter, “processing lag time”) can be analyzed equally, we performed RNA-Seq of RNA isolated from CHO-S cells (at 3, 5, 6, and 7 days of cultivation) with RNA protected at various times (15, 30, 45, and 60 min) following cell harvest, as shown in Fig. 1. High-quality total RNA was obtained from all samples (RIN values > 9.1). We prepared cDNA libraries from total RNA and obtained 11 to 42 million sequence reads per sample. The mean base quality score in the Phred scale was 35.68 indicating that the base call accuracy was above 99.97%. Sequence mapping rates to reference sequences ranged from 78.19% to 84.29%.

Shannon’s information entropy of RNA-Seq data of bach CHO-S cell culture was shown in Fig. 3. There were no significant differences in the value of Shannon’s information entropy between samples from the same cultivation day with different processing lag times. On the other hand, the value of Shannon’s information entropy decreased with cultivation. Decreasing Shannon’s information entropy during the cultivation is also seen in other organisms used in industry (unpublished data).

PCA (Fig. 4) and HCA (Fig. 5) were conducted using RNA-Seq data mapped onto the CHO-K1 RefSeq assembly (GCF_000223135.1) and CHO-K1 mitochondrial DNA (GCF_000055695.1). PC1 (cumulative contribution ratio (%) = 0.37) appeared to be associated with cultivation day (Fig. 4A), but neither PC1 nor PC2 (cumulative contribution ratio (%) = 0.10) were associated with the processing lag time (Figs. 4A–4E). The HCA dendrogram showed that samples of the same cultivation day clustered together; there were four clusters associated with cultivation day (Fig. 5). These results show that the gene expression was strongly associated with cultivation day, but the correlation between the gene expression and the processing lag time was not found. Even if processing lag times are different, global trends of gene expression were similar between samples from the same cultivation day.

Differentially expressed gene (DEG) analysis was conducted to explore the correlation between the processing lag time (15, 30, 45, and 60 min) and the gene expression. First, we compared the gene expression of RNA-Seq data between 15 min and each other time period (30, 45, and 60 min) on each cultivation day. The numbers of common DEGs among the four different days were 4 (15 min vs. 30 min), 10 (15 min vs. 45 min), and 23 (15 min vs. 60 min) (Fig. 6). The number of common DEGs increased as the processing lag time increased (Fig. 6, Table 1). All 23 common DEGs, with the exception of LOC103159497, were upregulated in 30, 45, and 60 min samples compared with 15 min samples (Table 1, Fig. S1). A total of 12 of the 23 common DEGs were orthologous to genes associated with apoptosis in the MGI database (Table 1).

PCA and HCA was then performed for genes identified as differentially expressed at least once in DEG analyses. In the PCA plot, PC1 and PC2 appeared to be associated with the processing lag time of each cultivation day (Fig. 7A). PCA for each cultivation day also appeared to be associated with the processing lag time (Figs. 7B–7E). In the dendrogram (Fig. 8), we identified three major clusters associated with cultivation day: cluster 1 included all day 3 samples; cluster 2 included all day 5 samples; cluster 3 included all days 6 and 7 samples. It was consistent with the PCA (Fig. 7) that samples of days 6 and 7 were not clearly separated in the dendrogram. Then, each sub-cluster was associated with the processing lag time. For example, cluster 3, which included days 6 and 7 samples, has three sub-clusters: cluster 3–1 included 15 and 30 min samples; cluster 3-–2 included 30 and 45 min samples; cluster 3–3 included 45 and 60 min samples. The influences of the processing lag time on the expression of transcriptome increased with cultivation day.

Discussion

To evaluate the influence of differences in the processing lag time of cell samples (i.e., time until the addition of RNA lysis buffer) on RNA-Seq, we examined gene expression changes when RNA lysis buffer was added at various times after the harvest of CHO cells. Regardless of the number of days of cell cultivation (3, 5, 6, and 7 days) or the processing lag time (15, 30, 45, and 60 min), we obtained high-quality total RNA and high-quality sequence reads from all samples. The results of the value of Shannon’s information entropy (Fig. 3), PCA (Fig. 4), and HCA (Fig. 5) showed that, for each cultivation day, the global trends of gene expression were similar between the samples with various processing lag times. In contrast, the expression levels of specific genes were affected by the processing lag time of the samples, despite the total RNA being of high quality. The number of common DEGs including apoptosis-associated genes increased as the processing lag time increased (Fig. 6, Table 1). In the dendrogram of the HCA performed for genes identified as differentially expressed at least once in DEG analyses (Fig. 8), each cluster consisted of several sub-clusters associated with the processing lag time. Especially, cluster 3 including the samples of days 6 and 7 separated into three sub-clusters associated with the processing lag time rather than cultivation day. Thus, we conclude that RNA lysis buffer should be added within 30 min after cell harvest. Otherwise, changes in the expression of some genes, especially genes associated with apoptosis, will occur.

Several studies have reported that the use of RNA samples with RIN values of around 5.0–7.0 does not negatively affect estimates of the gene expression profile (e.g., Fan et al., 2016; Shen et al., 2018) or de novo assembly (Kono et al., 2016). However, our results indicated that even if using high-quality RNA (RIN > 9.1) for RNA-Seq, the gene expression levels of specific genes were significantly different between samples that share the same cultivation day but have different processing lag times. For instance, three genes of the AP-1 transcription factor complex (Fos, Jun, and Atf4) were detected as DEGs (Table 1), with higher expression as the processing lag time increased. AP-1 transcription factors regulate the expression of target genes involved in a wide range of cellular processes, such as proliferation, survival, and death. Fos is regarded as an early response gene to both biochemical and mechanical stress in several cell types (Peake et al., 2000; Kim et al., 2017). In the case of CHO cells, shear stress (Ranjan et al., 1996), virus infection (Zachos, Clements & Conner, 1999), and cobalt chloride stimulation (Gong et al., 2001; Zou et al., 2001) have been known to induce AP-1 activation. Here, although the time spent on centrifuging and pipetting was equal for all samples (Fig. 1), the time placed in media or 1 × PBS at room temperature was longer for samples with the longer processing lag time, suggesting that expression of AP-1 transcription factors increased in response to 1 × PBS stimulation and/or the conditions (temperature and CO₂ concentration) outside the incubator.

Because the above results demonstrate that the processing lag time influences the expression level of specific genes, it is necessary to know not only the quality of the data but also how the experiment was performed when acquiring RNA-Seq data from SRA and investigating the expression level of a specific gene. However, detailed information on experimental operations, such as the processing lag time, is usually not described in manuscripts.

In recent years, there has been growing focus on the importance of reporting detailed methods, resulting in several protocol repositories have been begun, such as Nature Protocol Exchange (https://protocolexchange.researchsquare.com/), Protocols.io (Teytelman et al., 2016), MethodsX (http://www.sciencedirect.com/science/journal/22150161), Journal of Visualized Experiments (JoVE, https://www.jove.com/), and others. These repositories are an open-access, free platform for sharing, and discussing protocols and can be used during research and manuscript submission to facilitate reproducibility. However, protocol repositories usually do not show the working time of the experiment, even if it is a time-sensitive experiment like RNA isolation. Detailed protocols of this study including the processing time can be found at Protocols.io (see Materials & Methods).

Insufficient metadata will reduce the value of sequencing experiments by reducing the reproducibility of the study and its reuse for other analyses (Stevens et al., 2020). There are metadata standards for RNA-Seq such as the MINISEQE (Minimum Information about a high-throughput Nucleotide Sequencing Experiment), the ENCODE (Encyclopedia of DNA Elements) standards, or the International Human Epigenome Consortium (IHEC) metadata model. By following those metadata standards, or the checklist system for different types of data implemented by the European Nucleotide Archive (ENA), experimental details of RNA-Seq, such as the length of the RNA population being selected and whether rRNA populations were removed before library prep, can be described without missing. Furthermore, the BioSamples database at EMBL-EBI does not restrict the types of attributes used to describe sample metadata (Courtot et al., 2018), allowing submitters to choose attributes freely. However, those metadata standards for RNA-Seq do not have attributes corresponding to the working time of RNA isolation. The metadata associated with this study including the processing time of RNA samples described according to the ENCODE standards can be obtained from SRA Run Selector (see Material & Methods).

Although efforts have been made regarding the description to improve the reproducibility of experiments, it is difficult to describe the time required for individual steps and the number of pipetting operations in protocols or metadata. In the future, it would be ideal to employ humanoid robots that can use pipettes, vortex mixers, incubators, refrigerators, and centrifuges for experiments. A robotic crowd (or cloud) biology laboratory (RCBL; Yachie, Robotic Biology Consortium & Natsume, 2017) makes it possible to implement materials and methods complementary and robustly (Ochiai et al., 2020). The development of a humanoid robotic platform for sample processing would contribute to accurate recording and improve reproducibility of experiment operations.

Conclusions

In this study, we showed that differences in the time lag from harvesting cells to protecting RNA affect the expression level of specific genes in RNA-seq experiments. As discussed above, the data submitter needs to provide detailed information of the experiment together with the sequence reads. The database user needs to be able to investigate the validity of the data including potential problems with how the experiment was performed when acquiring HTS data from the database. Expression levels of the 23 DEGs identified in this study could potentially be used to help check whether the RNA of biological samples was properly treated.

Supplemental Information