The developmental effects of pentachlorophenol on zebrafish embryos during segmentation: A systematic view

Pentachlorophenol (PCP) is a typical toxicant and prevailing pollutant whose toxicity has been broadly investigated. However, previous studies did not specifically investigate the underlying mechanisms of its developmental toxicity. Here, we chose zebrafish embryos as the model, exposed them to 2 different concentrations of PCP, and sequenced their entire transcriptomes at 10 and 24 hours post-fertilization (hpf). The sequencing analysis revealed that high concentrations of PCP elicited systematic responses at both time points. By combining the enrichment terms with single genes, the results were further analyzed using three categories: metabolism, transporters, and organogenesis. Hyperactive glycolysis was the most outstanding feature of the transcriptome at 10 hpf. The entire system seemed to be hypoxic, although hypoxia-inducible factor-1α (HIF1α) may have been suppressed by the upregulation of prolyl hydroxylase domain enzymes (PHDs). At 24 hpf, PCP primarily affected somitogenesis and lens formation probably resulting from the disruption of embryonic body plan at earlier stages. The proposed underlying toxicological mechanism of PCP was based on the crosstalk between each clue. Our study attempted to describe the developmental toxicity of environmental pollutants from a systematic view. Meanwhile, some features of gene expression profiling could serve as markers of human health or ecological risk.

within a narrow range of origins. To accomplish this purpose, RNA sequencing is required to provide global information. This approach might simplify the abnormal organismal scenarios caused by the millions of toxic chemicals in environment and might be particularly beneficial for risk assessment and management.
In this study, we analyzed the transcriptomic effects of PCP treatments on zebrafish embryos using an Illumina HiSeq 2000 sequencing platform at 10 and 24 hours post-fertilization (hpf), representing the beginning and end of the zebrafish segmentation stage, respectively. The results suggested that PCP had significant impacts on multiple biological events in the developing embryos. Three different processes, metabolism, transmembrane transport and organogenesis, were chosen to ensure the integrity of causation. A mechanical framework that connected these events was proposed in an attempt to understand the occurrence and evolution of PCP's developmental toxicity from a systematic view.

Basic sequencing information.
Our exposure experiments were illustrated as Fig. 1. After exposure, the sequencing data produced approximately 57-72 M total reads from 6 samples, and their mapping/unique mapping rates were all above 90%. The original RNA-Seq data are available in the Gene Expression Omnibus database (GSE75921). On most chromosomes, the distribution of the detected transcripts would not have suffered from PCP exposure, and only chromosomes 3, 6, and 19 varied in some extent (Fig. S1).
The results of the differential expression analysis of the 21,227 detected transcripts are shown in Fig. 2a. A higher concentration and longer duration of exposure to PCP resulted in more differentially expressed genes, particularly those that were upregulated. Compared with the solvent control, 6,087 transcripts were affected in the 24H group (200 μ g/l PCP at 24 hpf), of which 3,593 were upregulated and 2,494 were downregulated. The 10L group (5 μ g/l PCP at 10 hpf) was the only treatment that included more upregulated genes than downregulated genes: 521 and 845 transcripts were up-and downregulated, respectively. The differentially expressed genes in the embryos treated with the two concentrations of PCP for two different durations are shown in a Venn diagram (Fig. 2b). A total of 5,895 genes were affected by only one of the four treatments: 3,663, 1,386, 493, and 353 genes, respectively. In comparison, only 55 genes were affected by all four treatments. We listed the top 10 upregulated mRNA transcripts in each of the 4 groups (Table 1) according to their fold changes (FC) and found that the genes in the two 10 hpf groups were completely different, while some of the genes in the two 24 hpf groups overlapped. Enrichment analysis based on the gene ontology (GO) definitions. Although the changes in gene expression seldom indicated an explicit biological impact on the organisms, a number of the enriched genes still  provided useful clues that should not be neglected. Here, we performed a functional enrichment analysis of the genes based on the GO definitions using three domains: biological process (BP), molecular function (MF) and cellular component (CC). The complete list of the significant GO terms is shown in Supplementary Table 1-4, as there were too many entries to list here.
To further filter the instructive terms from the abundant, affected terms, the GO trees of high-dose groups were graphed according to the hierarchical relationships of the terms, and only direct affiliations were accepted as valid (Fig. 3). If the toxic effects were primary and original, then we would expect their occurrence and development to be observed from a very early time point; therefore, they must continue to propagate and be the communal effects of the 10H and 24H groups. The BP terms tended to reflect PCP's effects on tissue and organ development, such as the kidney, eye, heart (and vessels), and central nervous system. The MF terms highlighted internal binding and transporter activity. One notable class was the binding activities to retinoids, including retinol and its core metabolites, retinal and retinoic acid. The filtered components of the CC terms included extracellular matrix, neuron projection, and mitochondrion in both groups and cell junction and proteasome in the 24 H group.

Enrichment analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway definitions.
The pathway analysis was an important addition to the GO analysis. Figure 4 describes the top 20 affected terms based on the definitions and annotations from the KEGG database, which were sorted by enrichment level. The characteristics of the KEGG definitions highlighted metabolism-related pathways, such as glycolysis and pyruvate metabolism. Moreover, some terms, such as ABC transporter, HIF1 signaling pathway, and ECM-receptor interaction, appeared frequently. The pie diagrams illustrate the relative proportions of the different types of terms (Fig. 5). The shares of the metabolic pathways varied from 7% to 38% in the four groups, corroborating a strong influence of PCP on embryonic metabolism. Among them, carbohydrate metabolism and its child term, glycolysis, were notable metabolic terms, which is consistent with the fact that many glycolytic genes exhibit tremendous fold changes. Many other downstream pathways involving lipids, amino acids, and other cofactors were also significantly regulated by exposure at 24 hpf likely resulting from alterations of the intermediates of carbohydrate metabolism. The digestive system terms, such as fat digestion and absorption and vitamin digestion and absorption, also had close connections with the metabolism terms.
A mechanical framework for the developmental toxicity of PCP. The proposed toxicological process of PCP is shown in Fig. 6. We assumed that all of the effects of pollutant exposure may have had a common origin and evolved gradually during development from a specific early stage (for example, the zygote stage). Because numerous reports have demonstrated that there are close and mutual interactions between the assorted biological events that were influenced in our experiments, some of the significantly altered processes were selected and connected based on their established interactions. Energy metabolism was regarded as the origin of the effects because their outcomes were closely related to the physiological state of the organism. Somitogenesis and eye formation were the major upregulated and downregulated features at 24 hpf, respectively, regardless of the PCP dosage. The transmembrane processes most likely acted as intermediaries between different events. See the Discussion section for a more detailed description.

Validation of the sequencing data.
Because only approximately 5.6% and 4.0% of differentially expressed genes in the 10 hpf and 24 hpf embryos, respectively, were affected by low doses of PCP, qRT-PCR was performed on embryos that were exposed to high doses of PCP to validate the sequencing data at these two time points. The selected genes are listed in Table 2 and include gapdhs, ldha (for glycolysis), atp5ia, cox7a2 (for oxidative phosphorylation), egln3 (for the hypoxia response), abcb11a, abcc8 (for ABC transporter), mespa, her1, cryba2b (for morphogenesis), kcnc3b, scn1ba (for ion channel), and gnb5b (for other functions). The definition of "significant" was only validated when p < 0.05 and |FC| > 2 were met simultaneously, thus aligning with the condition determined by sequencing. Using PCR amplification and quantification, the results confirmed that these genes had expression levels similar to those observed in the RNA sequencing.

Determination of PCP concentration.
A previous study showed that the actual toxicant concentrations decreased sharply upon exposure in a polystyrene microplate 10 . Hence, the actual nominal concentration of an initial 200 μ g/l PCP dose was determined at the beginning and end of the exposure period using chromatography ( Table 3). The measured concentration at the end of the exposure was 112.57 ± 2.95 μ g/l, which was slightly higher than the half-nominal concentration. This actual value could be observed in extreme conditions in the actual environment after degradation/absorption 2 , indicating that our study had some ecological significance.

Discussion
The affected processes used to construct a mechanical framework should meet the following requirements: (1) they are supported by some relevant and statistically significant terms at both time points and (2) the involved significant genes or gene families are necessary for the biological process. Based on these criteria, metabolism, transmembrane transporter, and organogenesis were the main components of the framework.
Affected energy metabolism-associated events. The role of energy is extremely fundamental to biological metabolism. Although OXPHOS provides approximately 90% of the energy, glycolysis has the advantage of generating energy rapidly and is efficient under some situations 11 . PCP exposure significantly affected glycolysis, OXPHOS, and their downstream carbohydrate, amino acid, and lipid metabolism pathways.
Glycolysis and hypoxia. The unusual activities of genes in the glycolysis and subsequent lactate fermentation pathways were the most remarkable characteristics of the zebrafish transcriptome, particularly in the 10H group, and were quite similar to those of the 8 hpf transcriptome in our previous study 5 . Both high-dose short-term and low-dose long-term exposure induced glycolysis activity but affected different genes. Most of the glycolytic genes, e.g., gapdhs, eno1a, pfkfb3, aldocb, and ldha, were highly active in the 10H group (Table S5) but were restored to their normal levels at 24 hpf. In contrast, some of the other glycolytic genes, e.g., aldob, pfkma, and gpia, were upregulated at 24 hpf. Two key genes that were treated as indicators of the Warburg effect, ldha and pkma, were only induced in the 10H group, suggesting that glycolysis may lead to potentially different consequences at the two stages.
Hypoxia and its core indicator, hypoxia-inducible factor-1α (HIF1α ), were highly correlated with the glycolytic process. Even in aerobic glycolysis, in which the oxygen supply was adequate, HIF1α expression fluctuated abnormally 12 . In addition to the glycolytic genes, most of the top 10 upregulated genes in the 10H group, such as egln3, egln2, slc2a3b, and igfbp1a, were thought to respond to cellular hypoxia. For example, igfbp1 overexpression was strongly induced by HIF1α and delayed zebrafish development 13 . Nevertheless, HIF1α expression itself remained fairly stable throughout the process, similar to that at 8 hpf 5 , and was even downregulated in the 24H group.
In light of the present knowledge, HIF1α may have been affected by the continuous, high expression of egln3 and egln2. These two transcripts encoded prolyl hydroxylase domain enzymes (PHDs) 3 and 2, respectively, which could be induced by hypoxia and then destabilize and degrade the HIF1α protein 14 . The restricted HIF1α expression may protect the zebrafish embryos from severe hypoxia in this experiment, although the actual roles of PHDs remain controversial. However, a dilemma of why the HIF1α -induced genes could be activated in the absence of HIF1α activation emerged, and it lacked a convincing explanation. OXPHOS. Traditionally, PCP has been recognized as a typical OXPHOS uncoupler. However, PCP exposure resulted in significant (24H) or mild (10H) inhibition of the embryonic OXPHOS, electron transport chain (ETC), and proton transport processes (Table S5). We noticed that in complex II, an alternative of complex I that  shares part of the ETC and the citric acid cycle, two hydrophobic subunits, sdhc and sdhd, which transferred the electrons to ubiquinone (coenzyme Q10), were downregulated after exposure, while two hydrophilic subunits, sdha and sdhb, which have no Q10 pool, were unaffected. The extensive therapeutic roles of Q10 in cancer, male infertility, heart failure, and other pathologies have been highlighted recently [15][16][17] , findings that may indicate underlying correlations between metabolism and chemical toxicity/diseases. OXPHOS often opposes glycolysis and hypoxia/HIF1α , particularly in cancer studies 18,19 ; however, in our experiments, they were not strongly correlated.
Retinoid metabolism. Similar to a "bagman", retinoids exhibited complicated relations with energy metabolism and the subsequent metabolic events. The metabolism, binding and signaling of retinoids were crucial, particularly for facial and cephalic morphogenesis, visual perception, and somitogenesis 20,21 . Retinoic acid, a crucial maternal signal, was primarily regulated by retinaldehyde dehydrogenases (RALDHs) and one of the cytochrome P450 isoforms, CYP26A1 20 . As shown in Table S6, PCP disrupted retinoid binding at both 10 and 24 hpf. The affected families were retinol-binding proteins, cellular retinoic acid-binding proteins, RALDHs (at 10 & 24 hpf), RA receptors (RARs), and retinoid X receptors (RXRs) (at 24 hpf). cyp26a1 upregulation (FC = 4.25 in the 24H group) echoed the alterations of the RARs/RXRs, suggesting that PCP has the potential to disrupt the entire retinoic acid signaling pathway.

Changes in the embryonic electrophysiological transport processes. Taking humans as an exam-
ple, approximately 10% of all genes are transporter related, which is consistent with their biological significance in maintaining the homeostasis of the organism 22 . Both the cargo they carried and the proteins themselves were useful for determining the underlying mechanism. In this test, the results of GO and KEGG pathway analyses highlighted the importance of the transmembrane transport processes in response to PCP exposure (Table S7).
Solute carrier family. Some members of the solute carrier (SLC) family were consistently altered from 10 to 24 hpf. Three members, slc4a1a, slc16a9a and LOC100150452, were the primary differentially expressed genes in the 10H group. The last two were upregulated and likely regulated carnitine and creatine transport 23,24 , both of which are critical metabolic agents that participate in ATP production and cellular energy homeostasis 25 . Slc4a1a encodes the Band 3 anion transport protein, the primary Cl-HCO 3 exchanger that regulates intracellular pH and chloride levels and that is downregulated in various cancer tissues 26 . Slc38a3 exhibited the highest expression levels among the dozens of altered genes in the 24H group, and its increased expression correlated with potassium restriction and metabolic acidification 27 . A conspicuous cluster in the 24H group included three transcripts from the SLC25 family, slc25a4, slc25a25a, and slc25a38a. The SLC25 family is known to be associated with mitochondrial and cytosolic metabolism, along with subsequent physiological processes, such as OXPHOS and gluconeogenesis 28 . Among them, SLC25A4 and SLC25A25, which transport ATP and phosphate, respectively, were significantly inhibited by PCP treatment. Most of the remaining genes were involved in amino acid transport.  Ion channels and synaptic transmission. The role of ion channels in early vertebrate development has been broadly described. However, other studies have already demonstrated the crucial functions of ion channels and pumps in early embryonic morphogenesis and hatching [29][30][31] . Here, we observed that PCP treatment influenced the genes encoding cation channels, including potassium, sodium, and calcium channels. Most of these channels were voltage-gated channels, with the exception of the partial potassium channels, which are "inwardly rectifying channels". Compared with the expression levels of the SLCs, the expression levels of the ion channels were relatively lower, except for cacna1fb, which was expressed at 10-fold higher levels than the other members. This gene participates in signal transmission in synaptogenesis and synaptic connectivity in the retina 32 and skeletal muscle 33 . In contrast, kcnk5 and kcnh2 are mainly expressed in the kidney and heart, respectively 34,35 . PCP also disrupted two members of the degenerin/epithelial Na + channel family, accn1 and accn5; the latter is a bile acid-sensitive subunit 36 .
Two key functions of ion channels are to assemble and to regulate neurotransmitter (NT) receptors. Therefore, as expected, our enrichment analysis obtained various NT-related terms, including glutamatergic, serotonergic, dopaminergic, cholinergic, and GABAergic synapses. In addition to ion channels, AMPA receptors (AMPARs), a type of ionotropic glutamate receptor, were affected by PCP exposure. Physiological and pharmacological studies have demonstrated that AMPARs participate in the development of mammalian and zebrafish central nervous systems 37 . These features also suggested that the neurobehavioral effect could be regarded as a marker for PCP's developmental toxicity.
ATP-binding cassette importers and exporters. The ATP-binding cassette (ABC) transporter family includes transmembrane proteins that pump exogenous toxicants and drugs out of the cell and that act as a barrier. In our test, many family members, including abcb11a, abcg2a, abcg8, abcc12, and the pathway term ABC transporters were the most significantly altered in the 24L group. Zebrafish lack the ABCB1/P-glycoprotein, the major xenobiotic efflux pump in human and rodents. Its physiological roles are executed by the homologs ABCB4 and ABCB5 38 , thus demonstrating that abcb5 was upregulated and reflecting the stress due to PCP exposure. However, the most impressive representation was from abcb11 (also known as BSEP/SPGP), which exhibited a 27.21-fold change in the 24H group. Although the existing research in human cell lines showed that ABCB11 was rarely involved in xenobiotic metabolism 39 , our result was consistent with the significant upregulation of the abcb11 gene after perfluorooctane sulfonate exposure 40 , suggesting that this exporter of endogenous substances has a certain role in exogenous exposure. Likewise, abcg2/BCRP, which is a high-capacity urate exporter 41 , was essential for the disposition and excretion of xenobiotics in rainbow trout 42 .
Effects on segmentation stages of zebrafish morphogenesis. Zebrafish embryos gradually become a complex organism after 14 h of segmentation stages mainly via the mesodermal differentiation and anteroposterior axis establishment. The rudiments of some primary organs form, including myotome/sclerotome (from somite), pronephros, neurula (which mature further into the brain, eye, and ear), and heart 6 . Although most of these tissues were influenced by PCP, two of the developmental events, somitogenesis and lens formation, were the most significant and, therefore, are discussed here.
Somitogenesis. Somitogenesis (the formation of somites from the mesoderm along the anteroposterior axis) is the most characteristic feature of the segmentation stage and was seriously affected by PCP (Table S8). The Mesp and Hairy protein families have central roles in this event 43 . Their members, including mespaa, mespab, her1, and msgn1, showed moderate changes at 10 hpf and vigorous upregulation at 24 hpf. The oscillations in the expression levels of these genes in the pre-somitic mesoderm, the so-called "Clock and Wavefront mechanism" 44,45 , were maintained during segmentation; however, in general, their expression levels exhibited a decreasing trend. Therefore, their upregulation implied a delay of developmental progress and likely affected the subsequent formation of bones and muscles. Because somitogenesis was accurately regulated by the competition between Wnt/ Fgf signals and RA signals, the evidence supporting this hypothesis was that wnt3a, fgf8a, while cyp26a1 (the gene encoding an enzyme to metabolize RA into bio-inactive form of retinoids) were upregulated.
Eye lens and other sensory organs. Zebrafish eyes were clearly recognizable in the embryonic head at 24 hpf, and the development of the visual system was likely considered a representative response to xenobiotics during early development. The most downregulated GO MF term in both 24 hpf groups was structural constituent of eye lens, composed of members encoding the crystalline-family proteins, which are the essential components of the eye lens. Crystallins could be classified into the α -, β -, and γ -subfamilies, and PCP primarily influenced the β -and γ -crystallins (Table S9). This phenomenon could be explained by the fact that α -crystallin accounts for less than 1% of the crystallins in larvae prior to 10 days post-fertilization (dpf) 46 . Meanwhile, the γ M-crystallins, which are essential for underwater vision, were the most abundant subunits 47 . The embryonic ocular lens is responsible for focusing on objects at different distances and protecting the eye from ultraviolet light; it enlarges from 24 to 72 hpf in zebrafish 48 . Therefore, damage to the lens could likely encumber subsequent eye development.
Even as early as 10 hpf, the enrichment analysis identified abnormally expressed genes in the eye, the ear, and the olfactory system. Photoreceptor cell morphogenesis was interrupted by BDE47 in 6 dpf zebrafish larvae in our previous study 49 . These facts reminded us that environmental pollution can cause remarkable damage to the eyes and other sensory organs of fish. The senses, particularly vision, guide most larval and adult behaviors; therefore, their impaired functions may have profound impacts on fish survival and reproduction in real ecosystems.
Systematic views of developmental toxicity. Explaining the massive data obtained from the high-throughput tests has become the greatest bottleneck in -omics fields. It is impossible for all of the existing enrichment terminologies to connect cross-category events, although inextricable links, such as glycolysis and hypoxia, have been identified between them. This reason might contribute to why the action mechanism of pollutants could not be clarified by enrichment analysis alone and why we proposed a hypothetical framework to understand the toxic mechanism.
A basic setting for the framework is that the biological events that are most heavily affected are also indicators for PCP exposure. From a systematic viewpoint, we tended to focus on embryonic toxicological phenomena that are derived from abnormal energy metabolism, because of their consistent performance during 8-10 hpf and because of their profound consequences. In addition to the ATP outcomes, OXPHOS influenced the cellular membrane potential, and glycolysis affected other downstream metabolic pathways through lactate/pyruvate and altered cellular acidity. These metabolic events also had inherent connections with cellular oxygen and hypoxia. The changes in glycolysis/carbohydrate metabolism have also been observed in large toxicological studies [50][51][52] but have rarely been investigated. The subsequent physiological conditions further influenced transmembrane transport. pH controls, enhances or inhibits many members of SLC and ABC transporter families, and ATP/ADP modulates the K currents from the voltage-gated potassium channels 53 . The metabolic products of glycolysis are known to be essential for subsequent metabolic pathways, such as lipid and amino acid metabolism.
At 24 hpf, the embryos have sufficiently matured to display morphological features. Thus, the toxic effects on organ formation could be the extrinsic and terminal parts of the framework. Moreover, the impacts of PCP on zebrafish organogenesis seemed to be associated with retinoid actions. Retinoids have two main active forms: all-trans retinoic acid and 11-cis retinal. Given the roles of retinoic acid signaling in somitogenesis (antagonist) and eye development (agonist) 54,55 , PCP was suggested to influence these two processes by repressing RA signaling. RA signaling has been shown to regulate the formation of the mesoderm and anteroposterior axis, and retinoid metabolism has been associated with other metabolic pathways. Previous studies have described possible interactions between retinoids and glycolysis/hypoxia, and other studies may be ongoing.
Here, we proposed a framework to study the mechanism of xenobiotic exposure from a systematic and progressive viewpoint. The framework connected disparate biological events, particularly physiological events (metabolism), to gene regulation events (organogenesis) to convert the disorganized transcription information into a readable mechanical story. If we attained our goal, we expected to identify common parts in the toxic mechanisms of millions of chemicals. However, our present study was challenging, and more details must be collected to fill in the gaps in this framework.

Methods
Zebrafish and toxicological exposure. Wild-type Tuebingen zebrafish (Danio rerio) were maintained in a recirculating filtration system using water treated by reverse osmosis (ESEN, China). The protocols for zebrafish feeding and spawning and water quality monitoring have been previously reported 49 . PCP (purity > 98%, Sigma, New Haven, CT) was dissolved in dimethyl sulfoxide (DMSO) as the vehicle (purity > 98%, Amresco, Solon, OH). The final concentrations of PCP were 0, 5, and 200 μ g/l, each of which contained 0.01% vehicle. The exposure started at approximately 1 hpf in a 28.5 °C incubator and was stopped at 10 and 24 hpf, respectively. Approximately sixty embryos in each treatment group were used for deep sequencing. For the quantitative PCR assay, the PCP treatments and DMSO controls were replicated 3 times. All animal protocols were carried out in accordance with guidelines for care and use of laboratory animals as approved by the Animal Ethics Committee of Tongji University.

Chemical determination.
To investigate the influence of the exposure duration on the PCP concentration, the actual concentrations of a nominal 200 μ g/l PCP solution were determined at the beginning and end of the exposure using an Agilent 1200 liquid chromatography system (Waldbronn, Germany) with a VW detector operating at a wavelength of 305 nm. The constituents of the mobile phase were methanol and 1 g/l ammonium acetate (V/V: 70/30), with a flow rate of 1 ml/min. The injection volume was 20 μ l. A TC-C18 reversed-phase column (150 mm × 4.6 mm, Agilent Technologies, Palo Alto, CA, USA) was used, and the column temperature was set to 25 °C. The separation was completed in 8 min.
RNA sequencing experiments, differential gene expression analysis and functional enrichment analysis. The experimental and analysis protocols have been previously reported 56 . The total RNA was extracted and purified using an RNeasy Mini Kit (Qiagen, Germany) according to the manufacturers' instructions. The RNA library was constructed, and sequencing was performed using an Illumina HiSeq 2000 platform (San Diego, CA, USA) after the RNA quality was examined using an Agilent TapeStation 2200 (Santa Clara, CA, USA). The clean reads from the raw sequencing data were mapped into the reference zebrafish genome Zv10 assembly using Tophat package v2.0, and further into human GRCH37.p13 genome by BLAST algorithm to achieve sufficient annotations. The mapped data were then normalized using the trimmed mean of the M-values (TMM) in edgeR package, and TMM counts over 0 were defined as detected. The statistical analysis for differential expression was performed using the R language tool of DEGseq, and "significance" was validated when p < 0.05 and |FC| > 2. To identify the statistical significance of the differentially expressed genes between the control and treatment groups, the gene functional enrichment analysis was then analyzed using Fisher's exact test and the definitions and annotations of the GO terminology and KEGG pathway databases. The significance threshold was set as p < 0.05, with the exception of the GO BP terms, which used p < 0.001. The GO trees were illustrated based on the hierarchies between the significant terms using Cytoscape 3.2.1, and only a direct relationship between the ancestors and children was adopted.
Quantitative real-time PCR and data analysis. After quality control of the extraction, the total RNA was reverse-transcribed into cDNAs with Superscript III Reverse Transcriptase (Invitrogen), random hexamers and oligo(dT) primers. SYBR Green I Quantitative PCR was performed on a 7500 Fast Real-Time PCR system (Applied Biosystems, Foster City, CA, USA). The PCR amplification mixture consisted of a 20-μ l final volume containing 0.2 μ l of Platinum Taq Polymerase, 2 μ l of PCR buffer (10×), 0.5 μ l of the forward primers, 0.5 μ l of the reverse primers, 1 μ l of SYBR Green dye (20×) and 1 μ l of the reverse transcription products. The PCR conditions were 40 cycles at 95 °C for 10 s, 60 °C for 30 s, and 70 °C for 30 s, and melting curves were generated after the last cycle. According to previous reports 57,58 and the transcriptomic information, elongation factor 1A (EF1A) was eventually chosen as the reference gene. The RT-PCR primers were designed using Primer Express software (Applied Biosystems) and are shown in Supporting Information Table S1. The threshold cycle (CT) values for each of the selected genes and EF1A were used to acquire equivalent amounts of RNA.