Comparison between the Amount of Environmental Change and the Amount of Transcriptome Change

Cells must coordinate adjustments in genome expression to accommodate changes in their environment. We hypothesized that the amount of transcriptome change is proportional to the amount of environmental change. To capture the effects of environmental changes on the transcriptome, we compared transcriptome diversities (defined as the Shannon entropy of frequency distribution) of silkworm fat-body tissues cultured with several concentrations of phenobarbital. Although there was no proportional relationship, we did identify a drug concentration “tipping point” between 0.25 and 1.0 mM. Cells cultured in media containing lower drug concentrations than the tipping point showed uniformly high transcriptome diversities, while those cultured at higher drug concentrations than the tipping point showed uniformly low transcriptome diversities. The plasticity of transcriptome diversity was corroborated by cultivations of fat bodies in MGM-450 insect medium without phenobarbital and in 0.25 mM phenobarbital-supplemented MGM-450 insect medium after previous cultivation (cultivation for 80 hours in MGM-450 insect medium without phenobarbital, followed by cultivation for 10 hours in 1.0 mM phenobarbital-supplemented MGM-450 insect medium). Interestingly, the transcriptome diversities of cells cultured in media containing 0.25 mM phenobarbital after previous cultivation (cultivation for 80 hours in MGM-450 insect medium without phenobarbital, followed by cultivation for 10 hours in 1.0 mM phenobarbital-supplemented MGM-450 insect medium) were different from cells cultured in media containing 0.25 mM phenobarbital after previous cultivation (cultivation for 80 hours in MGM-450 insect medium without phenobarbital). This hysteretic phenomenon of transcriptome diversities indicates multi-stability of the genome expression system. Cellular memories were recorded in genome expression networks as in DNA/histone modifications.


Introduction
When environmental conditions change abruptly, living cells must coordinate adjustments in their genome expression to accommodate the changing environment [1]. It is possible that the degree of change in the environment affects the degree of change in the gene expression pattern. However, it is difficult to completely understand the amount of change that occurs in the transcriptome, given that this would involve thousands of gene expression measurements.
A recent study defined transcriptome diversity as the Shannon entropy of its frequency distribution, and made it possible to express the transcriptome as a single value [1]. Dimensionality Reduction methods e.g. Principle component analysis and t-SNE had been used to transcriptome analyses [2,3], but biological meanings of value of results of these methods was unclear. Transcriptome diversity also reduces the dimensions of transcriptome data and the biological meaning of transcriptome diversity is clear. The first research on transcriptome diversity was performed with human tissues [1], and later research compared cancer cells with normal cells [4]. In plant research, transcriptome diversity has been used to compare several wounded leaves [5]. Our previous study indicated that silkworm fat-body tissues cultured for 90 hours had higher transcriptome diversity than intact tissues and that transcriptome diversity measured degree of cells development [6]. However, these studies considered only qualitative environmental changes. Here, we present comparisons of transcriptome diversity between cells cultured in vitro in media supplemented with several concentrations of phenobarbital, to investigate how the amount of environmental change (in terms of drug concentration) affects the amount of transcriptome change.

Effect of phenobarbital on transcriptome diversity
To investigate the effects of phenobarbital concentration on transcriptome diversity, we sequenced 15 transcriptomes from larval fat-body tissues exposed to phenobarbital. Freshly isolated tissues were cultured for 80 hours in MGM-450 insect medium, and then cultured for 10 hours in medium supplemented with 0, 0.25, 1.0, 2.5, and 12.5 mM phenobarbital. We measured the diversity of those transcriptomes.
Transcriptome diversity were 9.96, 10.39, 10.55 with 0 mM phenobarbital, 10.18, 10.36, 10.43 with 0.25 mM phenobarbital, 8.17, 7.57, 8.00 with 1.0 mM phenobarbital, 7.68, 7.33, 8.06 with 2.5 mM phenobarbital and 8.00, 7.52, 7.93 with 12.5 mM phenobarbital. Correlation between transcriptome diversity and phenobarbital concentration was -0.53 (t = -2.28, df = 13, p-value = 0.04, 95 percent confidence interval, -0.94 to -0.54, Pearson's product-moment correlation). Transcriptome diversity changed only between phenobarbital concentrations of 0.25 and 1.0 mM (Fig 1). Our hypothesis was that the amount of transcriptome change is proportional to the amount of environmental change (i.e., drug concentration). These results showed that there is no proportional relationship. In this research, we studied the effects of the entire range of phenobarbital concentrations that would be tolerated by fat-body cells. However, only two values of transcriptome diversity were recorded.

Transcriptome diversity responds to cis-permethrin
There is a possibility that the changing transcriptome diversity a phenobarbital-specific phenomenon. To discuss this, we cultured fat bodies in media supplemented with two concentrations of cis-permethrin, an insecticide that is hydrophobically opposite to phenobarbital. We sequenced two transcriptomes of fat bodies cultured in media supplemented with 0.25 and 2.5 mM cis-permethrin. Transcriptome diversities were calculated as 10.3 (0.25 mM cis-permethrin) and 8 (2.5 mM cis-permethrin) (S3 Fig). This result matches perfectly with those from the phenobarbital experiments. Therefore, the increase in storage-protein gene expression cannot be explained as a phenobarbital-treatment-specific phenomenon.

Hysteretic phenomena of transcriptome diversity
In this study, there is no proportional relationship between drug concentration and transcriptome diversity. This relationship looks non-liner, but it was not confirmed. To discuss that, we examined hysteresis between drug concentration and transcriptome diversity. Observation of the hysteretic phenomenon would indicate the existence of multi-stability. We cultured fat bodies in media supplemented with lower concentrations of phenobarbital after previous cultivation in media supplemented with a higher concentration of phenobarbital. It is becoming increasingly clear that many biological systems are governed by highly non-linear bi-or multistable processes, which may switch between discrete states, induce oscillatory behavior, or define their dynamics based on a functional relationship with the memory of input stimuli [7]. Hysteresis in molecular biology is known in a synthetic mammalian gene network [8]. There is a possibility that hysteresis is hidden in the non-proportional relationship between external stimuli and the transcriptome as a huge complex of gene networks. Topological changes are known in gene regulatory networks [9]. In response to diverse stimuli, transcription factors alter their interactions to varying degrees, thereby rewiring the network. It is possible that our results indicate hysteresis and multi-stability of the transcriptome.
To investigate hysteretic phenomena in transcriptome diversities, we sequenced six transcriptomes of fat bodies cultured in media supplemented with 0 and 0.25 mM phenobarbital for 10 hours, after 90 hours' previous cultivation (cultivation for 80 hours in phenobarbitalnon-supplemented MGM-450 insect medium, followed by cultivation for 10 hours in 1.0 mM phenobarbital-supplemented MGM-450 insect medium). We measured the diversity of those transcriptomes. Transcriptome diversities were 10.13, 9.75, 10,26 in 0 mM phenobarbital after the 1.0 mM phenobarbital experimental group, and 9.350, 10.22, 9.11 in mM phenobarbital after the 1.0 mM phenobarbital experimental group (Fig 1). Transcriptome diversity was determined not only by the drug concentration at that time, but also by previous drug concentrations. This is a hysteretic phenomenon, and a hysteretic phenomenon provides evidence of a bi-stable system. These results indicate multi-stability of the genome expression system.

Transcriptome diversity and differentially expressed genes
Determination of the drug concentration is critically important for in-vitro drug-exposure testing in preclinical toxicology. High drug concentrations can induce a radical transcriptome response, while mild concentrations can make it difficult to determine cell responses. Ideally, a drug concentration should be high enough to induce the desired main effect, while not eliciting too many side effects. If we performed differentially expressed gene analysis using several drug concentrations, we would obtain more than 100 differentially expressed genes from each comparison and several thousands of differentially expressed genes in total. We hypothesized that it would be possible to determine a perfect drug concentration by referencing transcriptome diversity as the index of the extent of transcriptome changes.
We compared transcriptomes with close transcriptome diversities and those with distant transcriptome diversities, focusing on differentially expressed genes. In this study, a comparison between the transcriptome of a control group (0 mM phenobarbital) and that of an experimental group treated with the lowest phenobarbital concentration (0.25 mM) represented the comparison of close transcriptome diversities. Indeed, only 29 differentially expressed genes were detected in the comparison between the control and the 0.25 mM phenobarbital-treated group. On the other hand, more than 1000 genes were detected in comparisons between the control group and the 1.0, 2.5, and 12.5 mM phenobarbital groups (1534, 2198, and 1302 differentially expressed genes, respectively). We further examined these differentially expressed genes.
It is known that insects have evolved mechanisms to detoxify, reduce their sensitivity to, or excrete insecticides [10]. Increased detoxification can occur by gene duplication of carboxylesterase, which cleaves the insecticide, or by transposon insertion, causing increased transcription of cytochrome P450, which hydroxylates the insecticide [10]. Point mutations in the target gene can reduce insecticide binding [10], while increased transporter activity leads to faster excretion from the cell [10]. Gene duplication of carboxylesterase and point mutations in the target gene have been detected using comparative genomics; however, increased transcription of cytochrome P450 and increased transporter activities can be detected using comparative transcriptomics. Phenobarbital is well known to induce additional cytochrome P450s and other drug-metabolizing enzymes [11]. The effects of hepatocyte growth factor on phenobarbital-mediated induction of ABCB1 [a member of the ATP-binding cassette (ABC) transporter superfamily] mRNA expression have been investigated [12]. Phenobarbital is a substrate for Pglycoprotein, which is an energy-dependent transmembrane efflux pump encoded by the ABCB1 gene [13]. Therefore, we focused on phenobarbital-related induction of cytochrome P450 and ABC transporter coding genes. In the current study, cytochrome P450 and ABC transporter genes were identified as differentially expressed genes in comparisons between the control and treatment groups. The various concentrations of phenobarbital induced two (0.25 mM phenobarbital), 12 (1.0 mM), 12 (2.5 mM), and 17 (12.5 mM) cytochrome P450 coding genes (Table 1). Marked elevations in cytochrome P450 expression levels were observed only in two cytochrome P450 genes (BGIBMGA001004 and BGIBMGA001005). These genes were detected from comparisons between the 0.25 and 12.5 mM phenobarbital groups. The ABC transporter coding gene (BGIBMGA007738) increased its expression only in the 0.25 mM phenobarbital experiment. In contrast, many ABC transporter genes in the !1.0 mM phenobarbital experimental groups decreased their expression levels. These results show that the comparison between transcriptomes with close diversities was useful in identifying genes induced by drugs.

Macroscopic similarity between transcriptomes
To discuss the biological origin of the decreasing transcriptome diversity, we compared the transcriptome diversities of cultured cells, cells cultured with phenobarbital, and intact cells. In a previous study, we compared the transcriptomes of intact and cultured silkworm fat bodies, and showed that fewer genes occupied more of the transcriptome in intact than in cultured fat bodies [6]. In intact fat bodies, storage-protein coding genes occupied more than half of the transcriptome and the diversity of that transcriptome was low (6.49). The transcriptomes of fat bodies exposed to higher concentrations of phenobarbital, with low transcriptome diversity, may have similarities to the transcriptome of the intact fat body.
To determine this, we sequenced two transcriptomes from intact larval fat-body tissue, and compared 18 transcriptomes using bar charts (Fig 2). The transcriptomes of intact larval fat bodies and those cultured in high concentrations of phenobarbital (>1.0 mM) were macroscopically similar. In these transcriptomes, storage-protein coding genes showed the highest expression levels, and these genes occupied more than one-third of their transcriptomes. Cultured tissue is thought to have initiated the process of abrogating its tissue-specific function by becoming independent from the donor and its identity as a part of an individual, living organism [6]. In this study, fat bodies exposed to high concentrations of phenobarbital recovered their tissue-specific character as part of an individual organism. It is conceivable that the culture including phenobarbital concentrations of !1.0 mM mimics the in vivo environment.

Discussion
The relationship between the degree of environmental change in response to external stimuli and the degree of transcriptome change was unclear. We compared the transcriptome diversities of silkworm fat-body tissues cultured with several concentrations of phenobarbital. We failed to find a proportional relationship, but cells cultured with phenobarbital concentrations of !1.0 mM had uniformly reduced transcriptome diversity. We also determined that the 500 genes with the highest relative frequency characterized transcriptome diversity by reducing the relative frequency of other genes (See Materials and Methods). Transcriptome diversity also demonstrated hysteresis. Multi-stability and hysteresis, phenomena widely known in nature [7] were also discovered in transcriptome responses. The concepts of gene expression state space and attractors had been introduced which provide a mathematical and molecular basis for an "epigenetic landscape" [14]. Discovering of multi-stability of genome expression countenance the concepts of gene expression state space and the concepts of "epigenetic landscape". Cellular memories were recorded in genome expression networks as in DNA/histone modifications. We used two drugs (phenobarbital and cis-permethrin) in this study. The decrease in transcriptome diversity in the high drug-concentration groups was observed for both agents. Storage-protein genes occupied more than one-third of all transcriptomes with low transcriptome diversity (<9). Storage-protein genes occupied more than one-third of all transcriptomes with low transcriptome diversity (<9). It could be argued that storage-protein coding genes added to expression following exposure to a highly concentrated drug. Phenobarbital and cis-permethrin are considerably different in hydrophilicity-phenobarbital is a hydrophilic molecule and cis-permethrin is hydrophobic-and it is interesting to consider that storage proteins might cope with the threat from these drugs in the same manner.
Comparative transcriptomics for in-vitro preclinical testing are widely performed as an alternative to animal tests. However, there is no quantitative method to determine the drug concentrations that should be used for in-vitro preclinical testing. Therefore, we are forced to determine appropriate drug concentrations by considering cell morphology and through various other examinations. Using transcriptome diversity, however, we can find the perfect drug concentration. It would be interesting to perform a follow-up study on changes in transcriptome diversity using human cells, since a culture-induced increase in transcriptome diversity has also been observed in a comparison between intact human liver tissue (Hij: 10.9) and HepG2 cells (Hij: 13.3), which is the cell line established from human liver carcinoma [15] and is widely used as an in-vitro model for liver cells [16].
For our drug-exposure experiments, we used cultured tissue with high transcriptome diversity. The transcriptome diversity of these tissues was decreased by exposure to higher drug concentrations; the decrease in transcriptome diversity induced by higher concentrated drug exposure was not the phenobarbital-specific response of cultured tissue. We hypothesize that cells have some control mechanisms with respect to transcriptome diversity, and that these mechanisms responded to the stimuli of higher drug concentrations. Such mechanisms might enable cells to work cooperatively in multicellular organisms.

Establishment of primary culture
The p50 strain of the silkworm, Bombyx mori, was grown on the fresh leaves of the mulberry, Morus bombycis. The female larvae of the fifth instar were aseptically dissected 3 days after the fourth ecdysis and the fat body was isolated. More than 100 chunks of tissue (approximately 2 mm 3 ) were excised from the fat bodies of 108 larvae. Those tissue particles were incubated in cell culture dishes (diameter, 35 mm; BD Biosciences, Franklin Lakes, NJ, USA) with MGM-450 insect medium [17] supplemented by 10% fetal bovine serum (BioWest, Nuaillé, France) with no gas change. The tissue was cultured for 80 hours at 25°C. The microbes were checked for infection using a microscope. Infection-free tissues were used in the following induction assays. No antibiotics were used in the assays to maintain the primary culture.

Chemicals
All of the chemicals used in this study were of analytical grade. Phenobarbital sodium (Wako Pure Chemical, Osaka, Japan) was dissolved in distilled water to make a stock solution, which was added to the medium to make final concentrations of 0.25, 1.0, 2.5, and 12.5 mM phenobarbital. cis-Permethrin (Wako Pure Chemical) was dissolved in acetone. This solution was diluted with three times its volume of ethanol just prior to mixing with the medium. The final concentrations of cis-permethrin were 0.25 and 2.5 mM.

Induction assay
The original medium was replaced with phenobarbital-or cis-permethrin-containing medium in the induction assays. The primary culture tissues were incubated with 0.25, 1.0, 2.5, or 12.5 mM phenobarbital or with 0.25 or 2.5 mM cis-permethrin for 10 hours. The primary culture tissues for hysteresis analysis were incubated with 0 or 0.25 mM phenobarbital for 10 hours after 90 hours' previous cultivation (80 hours in phenobarbital-non-supplemented MGM-450 insect medium, followed by 10 hours in 1.0 mM phenobarbital-supplemented MGM-450 insect medium). Induction assays were terminated by soaking the tissues in TRIzol LS reagent (Invitrogen, Carlsbad, CA, USA), and the tissues were kept at -80°C until analysis.

RNA isolation
Total RNA was extracted using TRIzol LS (Invitrogen) and the RNeasy Lipid Tissue Mini Kit (Qiagen, Hilden, Germany) following the manufacturers' instructions. Silkworm fat bodies (30 mg) soaked in 300 μL TRIzol LS were homogenized. The homogenates were incubated for 5 minutes at 25°C, and the samples were subjected to centrifugation at 12,000 × g for 10 minutes at 5°C. The supernatant was then transferred to a new tube and incubated at 25°C for 5 minutes. Chloroform (60 μL) was added to each sample, and the homogenates were shaken vigorously for 15 seconds and then incubated at 25°C for 2 minutes. Samples were subjected to centrifugation at 12,000 × g for 15 minutes at 4°C; subsequently, the aqueous phase of each sample, which contained the RNA, was placed into a new tube. An equal volume (120 μL) of 70% ethanol was added and mixed well. The samples were transferred to an RNeasy Mini spin column and placed in 2 mL collection tubes. The lid was closed gently, and the samples were centrifuged for 15 seconds at 8,000 × g at 25°C. The flow-through was discarded. A total of 700 μL RW1 buffer was added to the RNeasy spin column. The lid was again closed gently, and the samples were centrifuged for 15 seconds at 8,000 × g. The flow-through was discarded. A total of 500 μL RPE buffer was added to the RNeasy spin column, and the samples were centrifuged for 15 seconds at 8,000 × g. The flowthrough was discarded. A total of 500 μL RPE buffer was added to the RNeasy spin column, and the samples were centrifuged for 2 minutes at 8,000 × g. The flow-through was discarded. The RNeasy spin column was placed in a new 2 mL collection tube, and the old collection tube with the flow-through was discarded. The samples were centrifuged at 13,000 × g for 1 min. The RNeasy spin column was placed in a new 1.5 mL collection tube. A total of 30 μL RNase-free water was added directly to the spin column membrane. To elute the RNA, the samples were centrifuged for 1 min at 8,000 × g. The eluate from the previous step was then added directly to the spin column membrane. To elute the RNA, the samples were centrifuged for 1 min at 8,000 × g. The integrity of rRNA in each sample was checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Library preparation and sequencing: RNA-seq
Sequencing was performed according to the TruSeq single-end RNA-sequencing protocols from Illumina for Solexa sequencing on a Genome Analyzer IIx with paired-end module (Illumina, San Diego, CA, USA). A total of 1 μg total RNA was used as the starting material for library construction, using the TruSeq RNA Sample Preparation Kit v2. This involved poly-A mRNA isolation, fragmentation, and cDNA synthesis before adapters were ligated to the products and amplified to produce a final cDNA library. Approximately 400 million clusters were generated by the TruSeq SR Cluster Kit v2 on the Illumina cBot Cluster Generation System, and 36-65 base pairs were sequenced using reagents from the TruSeq SBS Kit v5 (all kits from Illumina). Short-read data have been deposited in the DNA Data Bank of Japan (DDBJ)'s Short Read Archive under project ID DRA002853.

Data analysis and programs
Sequence read quality was controlled using FastQC program (http://www.bioinformatics. bbsrc.ac.uk/projects/fastqc). Short-read sequences were mapped to an annotated silkworm transcript sequence obtained from KAIKObase (http://sgp.dna.affrc.go.jp) using the Bowtie program. A maximum of two mapping errors were allowed for each alignment. Genome-wide transcript profiles were compared between samples. All the statistical analyses were performed using R software version 2.13.0 and the DESeq package. The homology search and local alignments were determined using Blast2Go [18]. Sequence data from Homo sapiens intact liver (SRA000299) [19] and the HepG2 cell line (SRA050501) were obtained from the DDBJ Sequence Read Archive (http://trace.ddbj.nig.ac.jp/dra/index.html). Short-read sequences were mapped to human cDNA obtained from RefSeq (http://www.ncbi.nlm.nih.gov/refseq).
We calculated transcriptome diversities using a script.

Genes characterizing transcriptome diversity
Changes in transcriptome diversity are based on changes in gene expression. However, it is not clear which genes contribute to such changes. To overcome these problems, we analyzed changes in transcriptome diversity. A previous study described the transcriptomes of each tissue as a set of relative frequencies, Pij, for the ith gene (i = 1, 2, . . ., g) in the jth tissue (j = 1, 2, . . ., t); and then quantified transcriptome diversity using an adaptation of Shannon's entropy formula: P ij log 2 ðP ij Þ H ij ¼ À P g i¼1 P ij log 2 ðP ij Þ Each term [Pij log2 (Pij)] represents a monotonic increase if Pij is larger than 0 and smaller than e -1 . The gene with the largest relative frequency makes the biggest contribution to transcriptome diversity but, in the current study, the largest relative frequency of genes could not be larger than e -1 . To compare Pij between samples, we sorted log2 (Pij) in order of Pij (S1 Fig). Transcriptomes with high diversity (0 and 0.25 mM phenobarbital experimental group) and low diversity (1.0, 2.5, and 12.5 mM phenobarbital experimental group) divided clearly. This result suggests that genes with extremely high relative frequency determine transcriptome diversity, reducing the relative frequency of other genes. To verify this, we estimated the diversity of transcriptomes that were lacking the top 10, 20, 30, 50, 100, 200, 300, and 500 most-expressed genes (S2 Fig). When we removed the top 500 genes, there was no difference in transcriptome diversity. These results show that the 500 genes with the highest relative frequency characterize transcriptome diversity by reducing the relative frequency of other genes.