Combinatorial expression of cell cycle regulators is more suitable for immortalization than oncogenic methods in dermal papilla cells.

The immortalized cell is an essential research tool that uses robust growth properties for the functional investigation of gene products. Immortalized mammalian cells have mainly been established using three methods: expression of simian vacuolating virus 40 T antigen (the SV40 method); human papilloma virus-derived oncoprotein E6/E7 (the E6/E7 method); or combinatorial expression of R24C mutant cyclin-dependent kinase 4, cyclin D1, and telomerase reverse transcriptase (the K4DT method). However, it is unclear as to which method is optimal for an in vitro model. Here, we compared the biological characteristics and genome-wide expression profiles of immortalized human dermal papilla cells generated by the SV40, E6/E7, or K4DT method. To our knowledge, this is the first study to comprehensively compare expression profiles to determine the optimal immortalization method for maintaining the original nature of the wild-type cells. These data would be valuable to scientists aiming to establish new immortalized cell lines.

Introduction

Cultured cells cannot proliferate infinitely due to shortening of the telomere repeat sequence. This limitation is known as “Heyflick's limit,” named for the scientist who discovered this phenomenon (Hayflick, 1965). The process of overcoming this limitation is called “immortalization,” and the immortalized condition allows us to continue the cell culture indefinitely. Based on the various types of senescence such as the replicative, premature, stress-induced, oncogene-induced, tumor-suppressor induced senescence, the cell proliferation is limited for the primary cells (Wang et al., 2019). The first immortalized cell line was the HeLa cell line, which was established from human cervical cancer cells (Jones et al., 1971). These immortalized cells became an essential tool for investigating gene function and play a critical role in the screening of candidate chemicals for pharmaceutical drugs. For the immortalization of cells, two methods are the most popular: expression of the oncogenic protein SV40 large T antigen (Tevethia, 1984) and expression of human papilloma virus-derived E6/E7 (Hawley-Nelson et al., 1989). Cellular immortalization is tightly connected to cancer development and cell cycle control since loss of function of negative regulators of cell cycle turnover is required for efficient and infinite cell growth.

SV40 large T antigen is an oncoprotein and was isolated from polyomavirus simian virus 40 (Tevethia and Ozer, 2015). The SV40 oncoprotein induces cellular immortalization by bypassing the negative regulation of retinoblastoma protein (pRB) and the suppression of the p53 tumor suppressor protein (Ahuja et al., 2005). Immortalization using SV40 is also classified as a transformation of mammalian cells from their dominant cell growth capability. Furthermore, SV40 protein is known to bind to various types of transcriptional regulators, such as p300 and CREB-binding protein (Cho et al., 2001).

The E6/E7 oncoprotein was isolated from human papillomavirus (HPV). HPV-derived E6/E7 is also known to induce cervical cancer in humans (Zur Hausen, 2002). The target protein of E6 is the p53 tumor suppressor protein (Münger and Howley, 2002). Expression of the E6 oncoprotein causes the degradation of p53 protein through ubiquitination (Scheffner et al., 1990). Expression of the E7 oncoprotein induces the inactivation of pRB and accelerates the transcriptional activity of the transcription factor E2F1, resulting in the elevation of cell proliferation (Münger et al., 1989). In the SV40 and E6/E7 methods, immortalization is achieved through inactivation of two major tumor suppressor pathways, pRB and p53. Because they cause inactivation of the tumor suppressor p53, which is characterized as the “guardian of the genome” (Lane, 1992), the expression of SV40 or E6/E7 oncoproteins sometimes causes genome instability, such as polyploid formation and chromosome gain or loss. Loss of function of the p53 gene causes genomic instability as a result of impairment of the accurate replication ability of the genome, and loss of function of p53 is frequently observed in around half of human malignant cancers (Royds and Iacopetta, 2006), indicating the critical role of p53 for tumor suppression.

In 2011, Shiomi et al. showed that combinatorial expression of R24C mutant CDK4, cyclin D1, and telomerase reverse transcriptase (TERT) results in the efficient immortalization of human-derived myogenic cells (Shiomi et al., 2011). Furthermore, Shiomi et al. showed that these immortalized myogenic cells retained the ability to differentiate into mature cells. Senescent cells halt cell proliferation due to the accumulation of the senescence protein p16 (Brenner et al., 1998). The accumulated cytoplasmic p16 protein binds to endogenous CDK4, which results in decreased kinase activity of the CDK4-cyclin D1 complex. The reduced kinase activity of the CDK4-cyclin D1 complex induces hypo-phosphorylation of the pRB protein, which results in a decrease in the protein level of the transcription factor E2F1, which is free from pRB protein binding. However, expression of R24C mutant CDK4 and cyclin D1 induces constitutively elevated CDK4-cyclinD1 complex activity since p16 protein cannot bind to mutant CDK4 due to this point mutation (Shiomi et al., 2011). The expression of R25C mutant CDK4 and cyclin D1 constitutively enhances cell proliferation in various types of human cells, such as myogenic cells, fibroblasts, corneal epithelial cells (Fukuda et al., 2019a), and dermal papilla cells (DPCs) (Fukuda et al., 2020). An important characteristic of immortalization using R25C mutant CDK4 and cyclin D1 is the intact condition of the p53 tumor suppressor gene, which enables cell expansion without genomic instability. Immortalization using R25C mutant CDK4, cyclin D1, and telomerase reverse transcriptase was named the K4DT method based on the genes used. Furthermore, it was been shown that the K4DT method works in a wide variety of animal cells, such as bovine, swine (Donai et al., 2014), prairie vole, monkey (Kuroda et al., 2015), midget buffalo, rat, and megabat-derived cells (Tani et al., 2019). This method can even be applied to sea turtle-derived cells, based on the highly conserved amino acid sequence and function of cell cycle regulators in animal evolution.

In this study, we comprehensively compare the biological characteristics of immortalized cells generated using three different methods (the SV40, E6/E7, and K4DT methods) on the same cell sources. The fundamental question is which immortalization method is best able to retain the original nature of the wild-type cells, as measured by the whole transcriptome. This study is the first report to describe the comprehensive expression profiles of different immortalization methods, which will contribute to understanding the characteristics of immortalized cells and could accelerate the application of established cells as an in vitro model for human diseases.

Results

Biological characteristics of immortalized human dermal papilla cells generated using the K4DT, SV40, and E6/E7 methods

We obtained primary human DPCs from PromoCell (Heidelberg, Germany) through a local distributor (Takara Bio, Shiga, Japan). After the expansion of the primary cells, we introduced either enhanced green fluorescent protein-expressing retrovirus, simian vacuolating virus 40 (SV40)-expressing retrovirus (Fukuda et al., 2012), or human papilloma virus-derived E6/E7-expressing retrovirus. As shown in Figure S1, although the efficiency of gene introduction was only around a few percent, by selecting the successfully introduced cells using antibiotics (G418 for SV40, hygromycin for E6/E7), we successfully established human DPCs that expressed SV40 or E6/E7. We previously reported the establishment of immortalized DPCs with combinatorial expression of R24C mutant cyclin-dependent kinase 4 (CDK4), cyclin D1, and telomerase reverse transcriptase (TERT) (Fukuda et al., 2020). Hereafter, we will refer to immortalized cells with the combinatorial expression of R24C mutant CDK4, cyclin D1, and TERT as K4DT. In Figure 1A, we show the morphology of wild-type DPCs and DPCs established using the SV40, E6/E7, and K4DT methods. Although E6/E7 and K4DT showed uniform cell size, there was a variation in cell size in SV40 cells. We first checked for the successful genomic insertion of the introduced gene with polymerase chain reaction (PCR) analysis. As shown in Figures S2 and S34, wild-type cells did not show any genomic insertion, while SV40, E6/E7, and K4DT cells showed the successful insertion of the expected combinations of expression cassettes. Furthermore, we detected the protein expression of CDK4, cyclin D1, SV40, E7, and tubulin with western blotting. As shown in Figures 1B and S35, while wild-type cells did not show intense signals of detected proteins, SV40, E6/E7, and K4DT cells showed protein expression of the corresponding introduced proteins. We next analyzed the cell cycles of the wild type, SV40, E6/E7, and K4DT cells. Although wild-type, K4DT, and E6/E7 cells did not show any abnormalities, SV40 cells exhibited a broader G2/M peak compared to wild-type cells, suggesting polyploid abnormalities in SV40 cells (Figure 2).

Figure 1

Cell morphology and the detection of introduced genes in wild-type and immortalized dermal papilla cells

(A) Cell morphology of wild-type, E6/E7, SV40, and K4DT immortalized DPCs at passage 3.

(B) Detection of introduced gene products in western blots of wild-type, K4DT, SV40, and E6/E7 cells.

Figure 2

Cell cycle analysis of wild-type, K4DT, SV40, and E6/E7 immortalized DPCs

A representative histogram of cells is shown.

Whole-gene transcriptome analysis of wild-type and immortalized DPCs

To comprehensively compare the expression pattern of the whole genome, we carried out RNA-Seq analysis using the Illumina Hiseq platform (100 bp paired-end). We show the workflow of the sequencing in Figure 3A. After the removal of the adapter sequences, we carried out mapping of the sequences. The amount of obtained sequence reads was at least 22 M reads, which is enough to describe the whole transcriptome (Figure 3B) (Fukuda et al., 2019b). To evaluate the reproducibility of the data, we carried out RNA-Seq with triplicate biological replications. The quality of the sequencing reads was determined using the TASTQC program (Figures S3–S6). Almost all the read data existed within the green area, indicating that the quality of the reads was reliable. We mapped the obtained sequence to an up-to-date human reference sequence (GRCh38) with the STAR program. The mapping ratio of reads was more than 95%, suggesting that our mapping strategy was appropriate for detecting the expression profile (Figure 3B). The mapping ratio to human genome sequence (GRh38) was 96.2% (wild type 1), 96.7% (wild type 2), 96.3% (wild type 3), 96.5% (K4DT-1), 96.1% (K4DT-2), 95.8% (K4DT-3), 95.1% (K4DT AR-1), 96.3% (K4DT AR-2), 96.4% (K4DT AR-3), 96.3% (SV40-1), 96.2% (SV40-2), 95.9% (SV40-3), 96.8% (E6E7-1), 96.5% (E6E7-2), 96.6% (E6E7-3). We, furthermore, mapped sequencing reads to total sequencing library of human cDNA (Homo_sapiens.GRCh38.cdna.all.fa) with bowtie2 program. The overall alignment rate of all reads to cDNA was below. 89.79% (wild type 1, 26838188 reads), 90.28% (wild type 2, 25388984 reads), 90.11% (wild type 3, 27158793 reads), 90.74% (K4DT-1, 23884949 reads), 90.74% (K4DT-2, 25895869 reads), 90.14% (K4DT-3, 26615408 reads), 88.95% (K4DT AR-1, 29396085 reads), 90.30% (K4DT AR-2, 26536855 reads) 88.04% (K4DT AR-3, 24322764 reads), 90.40% (SV40-1, 23171547 reads), 90.49% (SV40-2, 31207937 reads), 90.71% (SV40-3, 39095454 reads), 92.24% (E6E7-1, 23119054 reads), 92.57% (E6E7-2, 26509457 reads), 91.38% (E6E7-3, 28433365 reads).

Figure 3

RNA-Seq analysis of wild-type, K4DT, SV40, and E6/E7 immortalized DPCs

For a Figure360 author presentation of this figure, see https://doi.org/10.1016/j.isci.2020.101929.

(A) Workflow of the RNA-Seq analysis.

(B) The number of sequence reads of each sample obtained after PRINSEQ (blue bars) and the mapping ratio (yellow bars) are shown as a bar graph.

(C) Dendrogram analysis of wild-type, K4DT, SV40, and E6/E7 immortalized DPCs.

(D) Three-dimensional principal component analysis (PCA) based on the expression profiling of wild-type, K4DT, E6/E7, and SV40 immortalized DPCs.

The complete list of expression counts for wild-type, K4DT, SV40, and E6/E7 DPCs is listed on Figshare (https://figshare.com/s/124ed9cd873719d6cb4d). We inputted the expression profile of the whole genome into TCC-GUI for downstream analysis. Due to the limitation of TCC-GUI, we used top 25,000 genes to draw the dendrogram analysis and principal component analysis (PCA) plot to explain the variance of the expression level. First, we carried out dendrogram analysis. As shown in Figures 3C and S7, the results of correlation plot and data matrix among all replicates showed three biological replicates showed reproducible results on the RNA-Seq detection (Figure S7). Furthermore, we processed the data with three-dimensional PCA. From wild-type cells, the distance to K4DT cells was the smallest, the distance to E6/E7 cells was moderate, and the distance to SV40 cells was the largest (Figure 3D). We also showed the detailed locations of three-dimensional PCA analysis with a movie file (Figure 360). Three-dimensional PCA analysis based on 25,000 genes showed cumulative proportion of the variance to be 65.5% (PC1, 33.26%, PC2, 18.964%, PC3, 13.315%), indicating that 3D PCA would be classified as the results of whole-gene analysis. From these data, we concluded that the gene expression profile of K4DT immortalized cells was closest to that of wild-type cells when compared to E6/E7 and SV40 immortalized cells.

Figure360.

An author presentation of Figure 3

Movie presentation of PCA analysis based on the expression profiling of wild type, K4DT, SV40, and E6/E7 cells.

Visualization of differentially expressed genes among wild-type and immortalized DPCs

Using TCC-GUI (Su et al., 2019), we analyzed the expression profiles of the 33,122 human genes (https://doi.org/10.6084/m9.figshare.12197691.v2).

During the analysis of differentially expressed (DE) genes, the edgeR program extracted 15,103 genes as the first candidates (p value of edgeR is less than 0.01). Before applying to the DAVID functional analysis, we carried out the selection of DE genes with different expression threshold (no filtering, more than 300, 500, 700, 900, 1000, 1500, 1700, 2000, 2500, 3000, 10,000 counts) under the p value less than 0.01. The number of genes selected by the above mentioned threshold is summarized in Table 1. The number of genes was changed from 15,103 genes with no filtering to 576 with at least 10,000 counts in any of the sample. The full list of no filtering genes (https://doi.org/10.6084/m9.figshare.13174148.v1), list of at least 300 counts (https://doi.org/10.6084/m9.figshare.13174142.v1), list of at least 500 counts (https://doi.org/10.6084/m9.figshare.13174151.v1), list of at least 700 counts (https://doi.org/10.6084/m9.figshare.13174154.v1), list of at least 900 counts (https://doi.org/10.6084/m9.figshare.13174157.v1), list of at least 1000 counts (https://doi.org/10.6084/m9.figshare.13174166.v2), list of least 1500 counts (https://doi.org/10.6084/m9.figshare.13174172.v1), list of least 1700 counts (https://doi.org/10.6084/m9.figshare.13174175.v1), list of least 2000 counts (https://doi.org/10.6084/m9.figshare.13174178.v1), list of least 2500 counts (https://doi.org/10.6084/m9.figshare.13174202.v1), list of least 3000 counts (https://doi.org/10.6084/m9.figshare.13174205.v1), list of least 10,000 counts (https://doi.org/10.6084/m9.figshare.13174229.v1) are given.

Table 1

Number of differentially expressed (DE) genes on different thresholds on wild-type, K4DT, SV40, E6E7 immortalized human dermal papilla cells

Filtering	Number of DE genes	File name
No filtering, p < 0.01 on EdgeR	15,103	f_list1_1.csv
300 counts on any sample, p < 0.01 on EdgeR	9836	f_list1_2.csv
500 counts on any sample, P < 0.01 on EdgeR	8651	f_list1_3.csv
700 counts on any sample, p < 0.01 on EdgeR	8036	f_list1_4.csv
900 counts on any sample, p < 0.01 on EdgeR	6750	f_list1_5.csv
1000 counts on any sample, p < 0.01 on EdgeR	6350	f_list1_6.csv
1500 counts on any sample, p < 0.01 on EdgeR	4829	f_list1_7.csv
1700 counts on any sample, p < 0.01 on EdgeR	4355	f_list1_8.csv
2000 counts on any sample, p < 0.01 on EdgeR	3786	f_list1_9.csv
2500 counts on any sample, p < 0.01 on EdgeR	3060	f_list1_10.csv
3000 counts on any sample, p < 0.01 on EdgeR	2524	f_list1_11.csv
10,000 counts on any sample, p < 0.01 on EdgeR	547	f_list1_12.csv

Based on the DE gene list with different expression thresholds (no filtering, more than 300, 500, 700, 900, 1000, 1500, 1700, 2000, 2500, 3000, 10,000 counts), we carried out the pathway analysis on each threshold. The results of top 5 significant pathways in each threshold are summarized in Table 2. The number of times listed within top 5 was highest in cell cycle (10 times), ribosome (9 times), pathway in cancer (5 times). From the viewpoint of association to chromosome instability, we selected cell cycle and pathway in cancer.

Table 2

List of top five DAVID pathway analyses on differentially expressed genes detected by different thresholds among wild-type, K4DT, SV40, E6E7 immortalized dermal papilla cells.

Filtering	Number of DE genes	DAVID ranking
		1st	2nd	3rd	4th	5th
No filtering, p < 0.01 on EdgeR	15,103	Pathway in cancer	Small cell lung cancer	Neurotrophin signaling pathway	Chronic myeloid leukemia	Endocytosis
300 counts on any sample, p < 0.01 on EdgeR	9836	Cell cycle	Colorecal cancer	Pancreatic cancer	Chronic myeloid leukemia	Pyrimidine metabolism
500 counts on any sample, p < 0.01 on EdgeR	8651	Cell cycle	small cell lung cancer	Spliceosome	Ribosome	Pancretic cancer
700 counts on any sample, p < 0.01 on EdgeR	8036	Ribosome	Cell cycle	Ubiquitin mediated proteolysis	Chronic myeloid leukemia	Spliceosome
900 counts on any sample, p < 0.01 on EdgeR	6750	Cell cycle	Ribosome	Spliceosome	Ubiquitin mediated proteolysis	Chronic myeloidleukemia
1000 counts on any sample, p < 0.01 on EdgeR	6350	Ribosome	Cell cycle	Spliceosome	Proteasome	Focal adhesion
1500 counts on any sample, p < 0.01 on EdgeR	4829	Ribosome	Cell cycle	Focal adhesion	Proteasome	Pathway in cancer
1700 counts on any sample, p < 0.01 on EdgeR	4355	Ribosome	Focal adhesion	Cell cycle	Proteasome	Pathway in cancer
2000 counts on any sample, p < 0.01 on EdgeR	3786	Ribosome	Focal adhesion	Proteasome	Cell cycle	Pathway in cancer
2500 counts on any sample, p < 0.01 on EdgeR	3060	Ribosome	Focal adhesion	Proteasome	Regulation of actine cytoskeleton	Pathway in cancer
3000 counts on any sample, p < 0.01 on EdgeR	2524	Ribosome	Focal adhesion	ECM-receptorInteraction	Regulation of actin cytoskeleton	Cell cycle
10,000 counts on any sample, p < 0.01 on EdgeR	547	Ribosome	Focal adhesion	ECM-receptor interaction	Pathogenic Escherichia coli infection	Arrhythmogenic right ventricular cardiomyopathy (ARVC)

ECM, Extracellular Matrix

We obtained a gene list that corresponded to the cell cycle in DAVID (112 genes). The expression levels of these 112 genes, after normalization, are presented as a heatmap in Figure 5A. As shown in Figure 4A, the normalized heatmap reveals that the difference between wild-type and K4DT cells is small, while SV40 and E6/E7 cells have major differences compared to wild-type cells. The raw expression counts of the 112 cell cycle-related genes in the KEGG of wild-type, K4DT, SV40, and E6/E7 cells are shown using box plots in Figures S8–S14. The threshold for DE genes was a change in expression count of 2X or greater, or 0.5X or less, in immortalized cells, and these DE genes were mapped into the KEGG cell cycle gene pathway map (Figure 5). As shown in Figure 6, although the affected genes in K4DT cells were restricted to cyclin D-CDK4 complex-related networks, SV40 and E6/E7 cells exhibited significant differences in various types of growth signals, such as c-Myc, cyclin E-CDK2, cyclin B-CDK1, and the APC/C-Cdc complex.

Figure 5

Summary of upregulated and downregulated genes among KEGG cell cycle-related genes

Arrowheads indicate upregulated or downregulated genes in K4DT, E6/E7, or SV40 cells. The genes with more than 200 counts and X2 increase or 1/2 decreased from the wild type were mapped.

Figure 4

Heapmap analysis

Heatmap of genes listed in cell cycle pathway (A) and pathway in cancer (B) in wild-type, K4DT, SV40, and E6/E7 immortalized DPCs. Red indicates high expression, and blue indicates low expression.

Figure 6

Summary of upregulated and downregulated genes among the KEGG pathway in cancer-related genes

Arrowheads indicate upregulated or downregulated genes in K4DT, E6/E7, or SV40 cells. The genes with more than 200 counts and X2 increase or 1/2 decreased from the wild type were mapped.

We also analyzed the pathway in cancer. The KEGG map of pathway in cancer contains 282 genes. The expression levels of these 282 genes, after normalization, are presented as a heatmap in Figure 4B. As shown in Figure 4B, the difference between wild-type and K4DT cells is relatively close when it is compared with that from SV40 or E6E7 cells. The raw expression counts of the 282 cell cycle-related genes in the KEGG of wild-type, K4DT, SV40, and E6/E7 cells are shown using box plots in Figures S15–S32. DE genes were mapped into the KEGG pathway in cancer map (Figure 6). Mapping on the pathway in cancer revealed that the c-Jun, c-Fos, c-Myc, which locate on the downstream of MAP kinase signaling, has unregulated in SV40 and E6E7 cells. Furthermore, expression levels of FGF and FGF-receptor were upregulated in SV40 and E6E7 cells.

Chromosome condition

Lastly, we carried out chromosome analysis of K4DT, SV40, and E6/E7 cells. We first detected the number of chromosomes in K4DT, SV40, and E6/E7 cells. As shown in Table S1, the K4DT cells had a normal human chromosome count of 46 in 50 mitotic cells. However, SV40 cells showed broad distribution of chromosome numbers, from 71 to 80, indicating that SV40 cells have major chromosomal abnormalities. Furthermore, the E6/E7 cells showed chromosome abnormalities in 10% of mitotic cells. We further analyzed the chromosome pattern of K4DT, SV40, and E6/E7 cells using the G banding method (Table S2). In K4DT cells, one mitotic cell showed translocation of chromosome 14 to 15, suggesting that the normal ratio of chromosome condition is 95%. However, SV40 cells had a normal ratio of 0%, and 20% of E6/E7 cells showed a gain or loss of chromosomes. From these results, we concluded that K4DT immortalization caused the fewest genomic abnormalities in DPCs. In contrast, the expression of SV40 caused major genomic instability, and the expression of E6/E7 induced minor or moderate genomic instability in DPCs (Figure 7B).

Figure 7

Chromosome analysis and summary of expression profiling of immortalized cells

(A) Chromosome analysis of K4DT, E6/E7, and SV40 cells using the G banding method. Small arrows indicate abnormalities detected by G banding.

(B) Summary of the expression profiling of K4DT, E6/E7, and SV40 immortalized DPCs. The distance between K4DT and wild-type cells is the shortest, that between E6/E7 and wild-type cells is moderate, and that between SV40 and wild-type cells is the greatest. Major genomic instability was observed in SV40 cells.

Discussion

In this study, we compared the whole genome expression profile among immortalized cells obtained by the expression of SV40, E6/E7 oncoprotein, or combinatorial expression of mutant CDK4, cyclin D1, and TERT (K4DT method). The establishment of immortalized cell lines from the same primary cell enabled us to compare the biological characteristics of the immortalization methods comprehensively. Based on the expression profiling of the whole genes, we concluded that K4DT is the most suitable method for maintaining the original nature of the primary cells. Based on previously established immortalized cells, such as HeLa or HEK293T cells, functional assays with these immortalized cells are essential to mimic in vivo. However, most scientists would wish to use immortalized cells with relatively intact characteristics that mimic the primary condition. Our manuscript first demonstrated that the combinatorial expression of mutant cyclin-dependent kinase 4 (CDK4), cyclin D1, and telomere reverse transcriptase is the most suitable method for immortalizing human dermal papilla cells while retaining their original nature, compared to the expression of SV40 antigen or human papilloma virus-derived oncoprotein E6/E7. Our data will be attractive to scientists trying to select the best method for in vitro cellular modeling of human and animal diseases. The anti-proliferative function of p16/tumor suppressor pRB has been listed as the first gate to cause replicative senescence in the cell replication (Wang et al., 2019).

In the pathway analysis with DAVID, cell cycle and pathway in cancer were selected for further analysis. However, the nucleotide excision repair pathway also plays a critical role on the maintenance of the genome. Based on this view point, we also mapped the DE genes on pathway (Figure S33). In good agreement with enhanced cell cycle turnover, upregulation of DNA polymerase complex including PCNA and DNA polymerase epsilon was detected. However, TFIIH complex-, XPE complex-, XPC complex-associated molecules did not show consistent change that suggest loss of function of repair complex. Potential abnormality in function of nucleotide repair will need to be analyzed about accuracy of genomic sequence of SV40 and E6E7 cells.

We also established K4DT immortalized DPCs expressing the androgen receptor (AR) gene using retrovirus gene transfer. Although the exact mechanism is unclear, the endogenous AR expression level became lower after several passages of DPCs (Kwack et al., 2008). In agreement with previous reports of AR suppression, we also detected downregulation of AR in DPCs immortalized using the K4DT method (Fukuda et al., 2020). We are now trying to analyze the whole expression data of K4DT immortalized cells expressing AR. The comparison between AR-expressing K4DTs and parent cells would contribute to our understanding of the androgen signaling pathway. A better understanding of the AR signal network would be beneficial for the development of anti-testosterone drugs to treat prostate cancer or androgenic alopecia, which is known to be affected by the activation status of AR.

Limitations of the study

In this study, we carried out the sampling from wild type around passage 3 to 4, while, the sampling from E6E7 cell was finished around passage 7 to 8. In K4DT cells, the sample collection was finished from passage 5 to 6 (see the section of method). Although we could not finish up the exact same number of passages for all cell lines, the accuracy of the passage number is extremely important for accurate evaluation. Furthermore, the accurate data for gene expression and biological triplicates might not be enough for the concise conclusion. Based on our data, more accurate and precise evaluation with larger biological duplication would be necessary to address the node genes that control the biological difference among immortalization method.

Resource availability

Lead contact

Further information and requests for resources and vectors should be directed to and will be fulfilled by the Lead Contact, Tomokazu Fukuda (tomofukuda009@gmail.com).

Data and code availability

The data sets with the next-generation sequencer that reported in this article are uploaded into the DNA Data Bank of Japan (DDBJ, and assigned BioProject accession SSUB015042. The output file for expression counts and lists of differentially expressed genes was listed in Figshare. The complete list of expression counts for wild-type, K4DT, SV40, and E6/E7 DPCs is listed on Figshare (https://figshare.com/s/124ed9cd873719d6cb4d). The full list of no filtering genes (https://doi.org/10.6084/m9.figshare.13174142.v1), list of at least 300 counts (https://doi.org/10.6084/m9.figshare.13174142.v1), list of at least 500 counts (https://doi.org/10.6084/m9.figshare.13174151.v1), list of at least 700 counts (https://doi.org/10.6084/m9.figshare.13174154.v1), list of at least 900 counts (https://doi.org/10.6084/m9.figshare.13174157.v1), list of at least 1000 counts (https://doi.org/10.6084/m9.figshare.13174166.v2), list of least 1500 counts (https://doi.org/10.6084/m9.figshare.13174172.v1), list of least 1700 counts (https://doi.org/10.6084/m9.figshare.13174175.v1), list of least 2000 counts (https://doi.org/10.6084/m9.figshare.13174178.v1), list of least 2500 counts (https://doi.org/10.6084/m9.figshare.13174202.v1), list of least 3000 counts (https://doi.org/10.6084/m9.figshare.13174205.v1), list of least 10,000 counts (https://doi.org/10.6084/m9.figshare.13174229.v1) are given.

Software and its versions were used for the data analysis.

FastQC, version 0.11.8, was used for quality check of raw FASTQ sequencing file. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/

PRINSEQ, version 0.20.4, was used for the removal of low quality reads. http://prinseq.sourceforge.net/

PEAT, version 1.2, was used for the removal of the adapter sequence. https://github.com/jhhung/PEAT.

STAR, version 2.6.1, was used for the mapping. https://github.com/alexdobin/STAR

featureCount, SUBREAD, release 1.6.5, was used for the expression counting. http://subread.sourceforge.net.

Bowtie2, version 2.3.4.3, was used for the mapping of reads to cDNA (Homo_sapiens.GRCh38.cdna.all.fa). http://bowtie-bio.sourceforge.net/bowtie2/index.shtml.

R package, version 4.0.3, was used for the downstream analysis. https://www.r-project.org.

TCC-GUI tool was used for the downstream analysis. https://github.com/swsoyee/TCC-GUI.

Materials availability

The established cell might be shared with the material transfer agreements of technology office and corresponding authors.

Methods

All methods can be found in the accompanying Transparent methods supplemental file.