text:
A systematic study on drug-response associated genes using baseline gene expressions of the Cancer Cell Line Encyclopedia

We have studied drug-response associated (DRA) gene expressions by applying a systems biology framework to the Cancer Cell Line Encyclopedia data. More than 4,000 genes are inferred to be DRA for at least one drug, while the number of DRA genes for each drug varies dramatically from almost 0 to 1,226. Functional enrichment analysis shows that the DRA genes are significantly enriched in genes associated with cell cycle and plasma membrane. Moreover, there might be two patterns of DRA genes between genders. There are significantly shared DRA genes between male and female for most drugs, while very little DRA genes tend to be shared between the two genders for a few drugs targeting sex-specific cancers (e.g., PD-0332991 for breast cancer and ovarian cancer). Our analyses also show substantial difference for DRA genes between young and old samples, suggesting the necessity of considering the age effects for personalized medicine in cancers. Lastly, differential module and key driver analyses confirm cell cycle related modules as top differential ones for drug sensitivity. The analyses also reveal the role of TSPO, TP53, and many other immune or cell cycle related genes as important key drivers for DRA network modules. These key drivers provide new drug targets to improve the sensitivity of cancer therapy.

Scientific RepoRts | 6:22811 | DOI: 10.1038/srep22811 gene mutation, copy number variation, and gene expression and built a drug-response prediction model using the selected biomarkers 8 .
However, cancer drug response mechanism is a very complex system that could be affected by many factors. Sex, in particular, can influence how the body handles a drug as well as the drug dose appropriate to the body 9 . Future clinical studies with patients using opioids for chronic pain should also include age as an important variable when assessing development of opioid tolerance 10 . Age is another important factor for the efficacy of drugs. It is known that there are more adverse drug reactions in the elderly than in the young, which might relate to the functional decline of clearing organs like kidney with age 11 . However, the specific genes and pathways involved in this process are not fully resolved.
Moreover, the DRA biomarkers may not function alone. Thus, it is of fundamental importance to identify not only individual gene markers, but also gene-gene interactions and modules associated with drug responses. For example, Chang et al. showed that pathway modules related to Ras-signaling and E2F transcription factors can be used to predict drug sensitivity 12 . By gene set enrichment analysis, it is also feasible to identify the key drivers or hub genes of a set of DRA genes in a regulatory or protein interaction network, the alternation of which will have substantial influence on the DRA gene set 13 . A number of network module methods have been developed and successfully applied in identifying the co-expression modules and key driver genes related to some diseases including Alzheimer's disease 13 , brain cancer 14 , and so on. However, to our best knowledge, this kind of analysis to drug sensitivity is till in its infancy.
In this paper, we developed a systems biology framework to identify gene expressions, co-expressions, and co-expression modules differentially changed with drug sensitivity. We then applied this framework to the Cancer Cell Line Encyclopedia (CCLE) 7 gene expression and drug response data and identified various DRA associated genes and modules. Moreover, we studied the effects of gender and age on DRA genes, and performed key driver analysis (KDA) 13 on differential functional modules (with drug sensitivity) defined by Gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways to identify key functional genes in the modules related to drug sensitivity.

Results
The CCLE data. The CCLE project is an effort to conduct a detailed genetic characterization of a large panel of human cancer cell lines 8 . It provides baseline gene expression profile of 20,069 genes for 504 human cancer cell lines collected from 24 tissue types and 21 cancer types. These cell lines are treated by 24 anti-cancer drugs including 17-AAG, AEW541, AZD0530, AZD6244, Erlotinib, Irinotecan, L-685458, LBW242, Lapatinib, Nilotinib, Nutlin-3, PD-0325901, PD-0332991, PF2341066, PHA-665752, PLX4720, Paclitaxel, Panobinostat, RAF265, Sorafenib, TAE684, TKI258, Topotecan, and ZD-6474. It is of note that very few cancer cell lines are present for some cancer types and tissues (e.g., bone and osteosarcoma), which might lead to biased conclusion. Thus, we removed cancer types and tissues with less than 5 cell lines. Finally, we obtained 323 cancer cell lines related to 17 tissues and 7 cancers. We plotted the distribution of sensitivity values for the 24 drugs on these cell lines in Fig. 1. For a cell line, the sensitivity of a drug is measured by the area under the dose-response curves (termed as activity area) 8 , which ranges from 0 to 7.8 in this study. In addition, we collected phenotype and annotation information of the patient from whom each cell line was generated, such as gender, age, and batch, to access their impact in drug sensitivity. The detailed information can be found in Supplementary Data S1.
Gene expression signatures of drug sensitivity. We first analyzed DRA gene expressions using a modified regression model (see Methods). Briefly, gender, age, tissue type, batch, cancer type, and top three principal components (PCs) for genotypes of cell lines are collected as potential confounding factors of drug sensitivity to cause gene expression variation. We listed the numbers of DRA genes (at false discovery rate (FDR) ≤ 0.1) for all 24 drugs in Table 1 and the detailed information for drug sensitivity association of all 20,069 genes in Supplementary Data S2.
As can be seen from Table 1, the number of DRA genes ranges from 2 to 1,226 (Lapatinib) across 24 drugs. Except for Lapatinib, there are also several hundred DRA genes identified for 17-AAG, AEW541, AZD6244, Erlotinib, Irinotecan, PD-0325901, Paclitaxel, Panobinostat, Sorafenib, TAE684, and Topotecan. Lapatinib is an orally active drug for breast cancer and other tumors 15 , and its treatment has been found to induce prevalent resistance in a couple of studies [16][17][18] . To further estimate how many false positives in our findings, we randomly permutated the sensitivity value of the cell lines for 1,000 times and counted how often we obtained more significant genes in permuted data than in the original one. In addition, we also counted the average number of significant genes for the 1,000 random permutations. The results were summarized in Table 1 under the column "Permutation". Except for AZD0530, L-685458, LBW242, Nilotinib, PF2341066, PHA-665752, PLX4720, and TKI258, the frequencies of obtaining more DRA genes in 1,000 permutation data than in original data for other drugs are less than or equal to 1. The results suggest that most of the identified genes are truly DRA. Due to the small numbers of DRA genes (≤ 20) in AZD0530, L-685458, LBW242, Nutlin-3, PF2341066, PHA-665752, PLX4720, and TKI258, and high number of false positives for Nilotinib and Sorafenib (i.e., 25 and 15.36), we removed them from further study unless stated.
As gene expression can be either positively or negatively associated with drug sensitivity and each could lead to different biological consequences, we divided DRA genes into positively and negatively associated genes (determined by the sign of the coefficient of the drug response term in the regression model). Except for Panobinostat, there are more positively associated DRA genes than negatively associated ones for other 13 drugs (Table 1). We visualized the expression pattern of DRA genes for 17-AAG in Fig. 2a and those for the other 13 drugs in Supplementary Fig. 1. As can be seen, all samples are divided into two large groups according to the hierarchical clustering (Euclidian distance with "Ward" measurement 19 ). The left side group consists of more samples with relatively larger sensitivity values, and is thus referred to as "sensitive group". The right group is referred to as "non-sensitive group". The clear separations of "sensitive" versus "non-sensitive" group are seen in all 14 drugs. We performed the Student's t-test between the sensitivity values of the two groups and the p-values are all smaller than 2.0 × 10 −7 (Table 2), indicating significant difference between the two groups. In addition, we plotted the expression levels of two sample DRA genes NQO1 and LOC344595 for 17-AAG (Fig. 2b,c), in which each dot represents a cell line. The Pearson correlation coefficients between drug sensitivity values and gene expressions for these two genes are 0.46 (p-value < 2.2 × 10 −16 ) and − 0.31 (p-value = 1.2 × 10 −8 ) respectively, which exhibits clear correlation between drug sensitivity and gene expression.
We next examined the overlap of DRA genes across 14 drugs. The results were summarized in Supplementary Data S3 and the top recurring genes were also plotted in Fig. 3a. Interestingly, many DRA genes occurred in multiple drugs. For example, DNAJB2 and ADCK were inferred to be DRA for 6 drugs. DNAJB2 (i.e., DnaJ (Hsp40) homolog, subfamily B, member 2) encodes a protein similar in sequence to bacterial DnaJ protein and the yeast homologs (http://www.ncbi.nlm.nih.gov/gene/3300). DnaJ proteins are differentially expressed in human tissues and demonstrate the capacity to function to both promote and suppress cancer development by acting as chaperones for tumor suppressors or oncoproteins 20 . Loss of expression for genes in this family confers resistance to chemotherapeutic agents used in treating ovarian cancer 21 . In addition, ADCK3 gene encodes a mitochondrial protein similar to yeast ABC1, which functions in an electron-transferring membrane protein complex in the respiratory chain. Expression of this gene is induced by the tumor suppressor p53 and in response to DNA damage, and inhibiting its expression partially suppresses p53-induced apoptosis 22 . These genes might be involved in some fundamental mechanisms of multi-drug targeting process.
DRA genes for drugs share common biological functions. The biological functions of positively and negatively associated DRA genes were annotated using David tools 23 . We plotted the word-cloud maps of the enrichment of DRA genes for two representative drugs Paclitaxel and Topotecan in Fig. 3b,c respectively and showed the enrichment of DRA genes for the rest drugs in Supplementary Data S4. The most enriched GO term for both drugs is membrane-enclosed lumen with FDR for Paclitaxel and Topotecan being 1.2 × 10 −49 and 5.8 × 10 −23 respectively. Membrane-enclosed lumen encompasses endoplasmic reticulum lumen, nuclear envelope lumen, and so on. A couple of studies have shown that changes in endoplasmic reticulum lumen environment affect cell sensitivity to apoptosis 24,25 and are thus important for cancer treatment 26 . Other top enriched terms include GO:0007049 cell cycle with FDR being 5.1 × 10 −16 and 2.6 × 10 −12 and GO:0003723 RNA binding with FDRs 1.6 × 10 −22 and 3.9 × 10 −11 respectively for the two drugs. One of the main mechanisms of anti-cancer drugs is to induce cell apoptosis 27 . In fact, there are numerous instances of cell-cycle mediated drug resistance 28 , however the underlying mechanisms are not fully known. The DRA genes related to cell cycle might shed some light on these mechanisms. In addition, the relationship between RNA-binding and drug sensitivity (or resistance) has also been widely studied 29,30 .
By an overview of all drugs, we identified a few functions significantly enriched in multi-drugs. For example, GO:0005886 plasma membrane, GO:0000087 M phase of mitotic cell cycle, GO:0000279 M phase, and GO:0000280 nuclear division are significantly enriched (at FDR ≤ 0.05) in four drugs, supporting their vital roles in drug response (Fig. 3d). The relationship between plasma membrane and drug sensitivity has been widely studied. For example, several cell membrane transporter proteins including multi-drug resistance protein 1 (MDR1; also known as P-glycoprotien and ABCB1), MDR-associated protein 1 (MRP1; also known as ABCG2), and  Table 1. Number of sensitivity-associated genes in 24 drugs. a Number of significant genes (at FDR ≤ 0.1) after adjusting gender, age, tissue, batch, cancer types, top three genotype principal components (PCs) as confounding factors. Columns "Positive", "Negative", and "Overall" list the number of positively-regulated, negatively-regulated and overall sensitivity-associated genes. b Information of sensitivity-associated genes (at FDR ≤ 0.1) in 1,000 permutation runs. Column "Frequency" lists the frequencies of identifying equal or more significant genes in the permuted datasets than those in the original one for 24 drugs; Column "Average Number" lists the average number of sensitivity-associated genes in 1,000 permutation runs. The blank in the table represents the drugs with the small numbers of DRA genes (≤ 20) and higher false positive rates judged by the permutation test.
Scientific RepoRts | 6:22811 | DOI: 10.1038/srep22811 breast cancer resistance protein (BCRP; also known as ABCG2) are known to cause resistance in common chemotherapeutics by promoting drug efflux 31 . Zaman et al. identified that the protein MRP, which is known to be associated to multi-drug resistance, is a plasma membrane drug-efflux pump 32 . Furthermore, as the hallmark of cancer revolves around cell-cycle deregulation, it is not surprising that antimitotic therapies are effective against the abnormal proliferation of transformed cells 33 . Besides the known functions, we also identified a few less known DRA functions, e.g., spindle. The results show that drug response is a very complicated process involving a lot of biological activities.
DRA gene patterns share certain similarity in male versus female. As the sample sizes of male and females are quite different (Supplementary Table 1), we first tested the effect of sample size on detecting DRA genes. Specifically, we randomly selected  10, 20, , 180 cell lines (each for 100 times) and calculated the average number of DRA genes for each sample size using the modified regression model (see Methods). The results for 9 drugs were shown in Fig. 4a and the full results for 14 drugs were summarized in Supplementary Data S5. As can be seen, the detection power varies significantly with sample size for most drugs. Thus, instead of directly comparing the number of significant DRA genes, we performed the Fisher's exact test 34 on overlapping genes in the top 100, 200, . . . , 5000 significant genes inferred from male and female samples separately. The results were shown in Supplementary Data S6. Interestingly, opposite patterns were identified for different drugs. For 11 (out of 14) drugs tested, the overlap becomes significant (at p-value ≤ 0.05) if more than 500 top genes are selected from male and female groups. The results suggest that male and female might share a lot commonality on the response In (a), each row indicates a gene and column indicates a sample; the heat-map colors represent gene expression with red for high expression and blue for low expression. We also added a side bar at the top to indicate sensitivity value with dark green for low values and yellow for high values. In (b,c), X-axis represents sensitivity and Y-axis represents gene expression level. Pearson-R value in the title represents the Pearson correlation coefficient between gene expression and sensitivity across all samples.
process of these drugs. However, the p-values are not constantly significant for the remaining 3 drugs. More differences (between male and female) were found for PD-0332991, RAF265, and ZD-6474. We performed an    Table 2. There is only one drug PD-0332991 targeted purely on gender specific cancers, i.e., breast and ovarian cancer. Interestingly, PD-0332991 is also one of the three drugs showing to have gender effect according to our regression model. The rest two drugs, i.e., RAF265 and ZD-6474, though not designed specifically for gender specific cancers, such as prostate or ovarian cancer, also have some sort of gender specific effect for part of their targeting cancers. For example, RAF265 is reported to inhibit the growth of human melanoma tumors 35 and men have much worse melanoma survival than woman 36 . ZD6474 is an inhibitor of VEGFR-2, which is designed for controlling the growth of lung metastasis and pleural effusions in human non-small cell lung cancer 37 . There are also clinical evidence showing that ZD6474 has potential role in the targeted therapy of breast cancer 38 . Men and women are shown to be different in both formation and prognosis of lung and breast cancers 39 . We also selected 9 drugs and plot their overlap in DRA genes between male and female in Fig. 4b. In summary, the two genders share common mechanisms for cancer drug responses for many drugs. However, there might be two patterns of DRA genes between genders.
Age effect of drug sensitivity. To investigate the impact of age in drug response mechanisms, we divided the samples evenly into two groups (namely young group and old group) by chronological age and repeated our analyses on these two groups. We listed the detailed information including sample sizes, numbers of significant DRA genes, and the age to separate the young and old for 14 drugs in Supplementary Table 3. As can be seen, the median age is 53 or 54 for all the drugs, which is not surprising since cancer is mainly an aging disease. The DRA genes show quite different patterns between young and old groups. Except for Irinotecan, RAD265, and ZD-6474, there are several to several hundred DRA genes for young group. However, the numbers of DRA genes for old group are close to zero for most drugs. Interestingly, there is only one significant DRA gene for 17-AAG in old group, however it is also identified as DRA for young group (with Fisher's exact test p-value for overlap being 2.4 × 10 −3 ). Similar scenario was observed for Topotecan (p-value 2.0 × 10 −3 ). It is known that age has a huge effect on the sensitivity of drugs and most drugs are more sensitive to young people 40 . Our results show that this difference might due to loss of function of some genes for the elder. However, we are fully aware that the small sample size should contribute at least in partial to the difference.

Drug-sensitivity associated gene modules and their key driver genes.
We have considered the relationship between individual gene expression and cell line-drug sensitivity. However, the sensitivity of a drug might relate to many genes and their interactions, thus it will be helpful to view the relationship from network perspective. Similar to Staunton et al. 41 , for each drug we divided the cell lines into three groups, namely sensitive group, indeterminate group, and resistance group. Specifically, cell lines with sensitivity values at least 0.8 standard deviations greater than the mean for a drug were defined as resistant to the drug; those with 0.8 standard deviations below the mean were defined as sensitive; and cell lines with sensitivity values within a window of 1.6 standard deviations around the mean were considered as indeterminate and were eliminated from analysis. To perform the differential network module analysis, we considered the biological modules defined by Gene Ontology 42 and Kyoto Encyclopedia of Genes and Genomes 43 (see Methods). It is of note that we also tried to define modules based on gene co-expression using WGCNA 14 . However, the modules obtained by sensitive and resistance group are not consistent and the sample sizes are also relatively small in each group. Thus, we adopted modules in GO and KEGG, which represent more general gene modules. Specifically, for each module, we first We then applied modular differential connectivity (MDC) analysis 13 to quantify co-expression difference between the resistant and sensitive states. MDC calculates the ratio between the average connectivity for gene pairs in resistant group and that for gene pairs in sensitive group. It is a continuous measure ranging from 0 to infinity. A module with MDC larger than (less than) 1 gains (loses) connectivity when changing from sensitive state to resistance state. The significance of MDC is estimated by a permutation study on cell lines (see Methods). We listed the MDC and significance for each GO and KEGG term in Supplementary Data S7. For a better view, we also showed top GO terms (with FDR 0 and extreme MDC values) for drug 17-AAG in Table 3. The top modules gaining and losing connectivity include GO:0032371 regulation of sterol transport, hsa04721 Synaptic vesicle cycle, and GO:0017156 calcium ion-dependent exocytosis. Cancer resistance proteins (e.g., ABCG and PCRP) have been implicated in the transport of sterols 44 , and exocytosis is one of the known mechanism related to multi-drug resistance in cancers 45 . In addition, modules related to cell membrane, drug metabolic process, synaptic transport, and immune systems are also in the top of the list.
For an overview of the MDCs of each GO term and KEGG pathway in multiple drugs, we listed their rank for each drug in Supplementary Tables 8 and 9 respectively and further ordered them according to the sum of the ranks across all drugs. Interestingly, we find that the GO term GO 0000003_reproduction ranks first at 7 out of 15 drugs tested. It is known that the power in inhibiting cell growth and reproduction is highly associated with cancer drug sensitivity 27 . In addition, consistent with the function analysis for DRA genes, modules related to cell cycle (e.g., cell cycle checkpoint, G1/S transition of mitotic cell cycle, and mitotic cell cycle) rank at the top.
For each significant module, it is very important to identify up-stream genes or hub genes directly interact with most genes in the module. The so called key driver genes or hub genes hold the key for reducing drug resistance and improving drug efficacy. They are usually identified by a few computational methods based on gene regulatory networks or protein interaction networks, e.g., key driver analysis (KDA) 13 . We selected a few significant modules for KDA analysis 13 . Specifically, we first constructed a protein-protein interaction (PPI) network using the interaction confidence level 0.4 in STRING database (version 10). We then mapped the genes in selected differential modules into the PPI network, and then retrieved the sub-network containing all the genes in the 2-neighborhoods of the mapped genes. By 2-neighborhoods, we mean the genes having a distance less than or equal to 2 with the mapped genes. We then identified the key driver genes (hub genes) of the mapped modules genes in the sub-network using KDA analysis 13 . It is worth noting that usually the key driver genes are inferred using a directed regulatory network like Bayesian network, however the key driver genes inferred from un-directed PPI network are also proved to be quite useful 46 .
We used module GO:0002444 myeloid leukocyte mediated immunity for drug Paclitaxel and module GO:1901992 positive regulation of mitotic cell cycle phase transition for drug Topotecan as two examples. The two networks were plotted in Fig. 5 using Cytoscape, in which the key drivers were in red and other genes in grey and the nodes with large size are more important. The top key driver genes for module GO:0002444 myeloid leukocyte mediated immunity are VEGFA and TSPO. VEGFA is a member of the PDGF/VEGF growth factor family, and has been identified to be related to cancer drug sensitivity in a couple of studies, e.g. 47 . Translocator protein (TSPO) is an 18 kDa high affinity cholesterol-and drug-binding protein found primarily in the outer mitochondrial membrane 48 . Interestingly, TSPO is also an important key driver gene for module GO:1901992 positive regulation of mitotic cell cycle phase transition. In addition, other genes including UBC, PCNA, and TP53 were inferred to be top key drivers. In fact, somatic mutations of TP53 gene are one of the most frequent alterations in human cancers. TP53 mutations are also potential prognostic and predictive markers, as well as targets for pharmacological intervention 49 .

Discussion
In this work, we presented a holistic view of the drug response related molecular markers including DRA genes and network modules for multiple drugs using the CCLE data. Both functional enrichment analysis on DRA genes and DRA network module analysis identified cell cycle as an important common function for multiple drug responses. It is known that many anticancer drugs functions by controlling the cell cycle which has been  Table 3. Top modules that gain or lose connectivity between sensitive and resistance groups for drug 17-AAG.
Scientific RepoRts | 6:22811 | DOI: 10.1038/srep22811 altered in human cancer 50 . In addition, our key driver analysis also infers TSPO, TP53, and many other immune or cell cycle related genes as key driver/hub genes of DRA modules. TSPO is important in binding small molecule drugs, cholesterol, and porphyrins 48 . The role of P53 in drug response has also been extensively studied 49,51 . For example, the p53 tumor suppressor are required to provoke cell apoptosis by DNA damage, which is crucial to the drugs inducing cancer cell apoptosis 51 . Destroying normal cell cycle, resisting apoptosis, and evading immune system are known hallmarks of cancer 52 . Interestingly, we also found the link between DRA genes and other cancer hallmarks. For example, GO:0000003 reproduction ranks first for 7 anticancer drugs among 2,846 GO terms (Supplementary Dataset S7), while the most significant hallmark of cancer is that the cancer cell can stimulate their own growth 52 . The modules related to other cancer hallmarks like genomic instability and inflammation are also rank as the top DRA modules. Thus the mechanisms of cancer drug sensitivity should be closely related to cancer hallmarks. In the future, it will also be interesting to study the potential of DRA genes or their co-expression patterns in predicting drug sensitivity by applying regression models. Our study revealed two different patterns of drug response between two genders. The DRA genes were significantly overlapped for drugs targeting cancers that are common in both male and female, suggesting that there might be some fundamental drug response mechanisms conserved. However, very different drug response genes were identified for drugs only targeting gender specific cancers (e.g., breast cancer). Female sex has been shown to be a risk factor for clinically relevant adverse drug reactions. The female specific DRA genes might be the key to understand why female usually need more drug dosage even though they have small weight in general. Similarly, age is a main factor for drug response and many studies have shown that people in different age differ in drug response 11 . Our study showed much more DRA genes for young people than for old people, which suggests that some DRA genes might lose their functions with the increase of age.
We are fully aware that there might be a few limitations of this work. For example, the sample size is not large enough in especially the age and gender studies, which contributes at least partially to the differences in DRA genes identified. In addition, though we removed many confounding factors to reduce the potential false-positive DRA genes, there could be other factors unmeasured. An alternative way is to apply surrogate variable analysis 53 to estimate these variables. However, this analysis usually estimate too many surrogate variables and it is not easy We use node shape to denote whether the node is in the module: (1) rectangle represents gene in module; (2) circle represents neighboring gene in the PPI network. We use color to denote whether the node is a key driver: (1) red represents key driver; (2) grey represents other gene. to determine a realistic number of variables. In addition, the estimated variables are not as reliable as the observed ones like gender and age. Thus, we prefer direct confounding factors as used in this study.
In addition, we used genes in GO and KEGG as modules to perform the differential module analysis due to the small sample sizes in sensitive and resistance groups. An alternative way is to define co-expression modules based on gene co-expression 13 , which is a more data-driven approach. It should also be better to construct Bayesian regulatory networks to infer the key driver genes, which considers the direction in regulation 13 . However, a large sample size is required for constructing robust co-expression and regulatory networks and we will try to integrate more samples in the future to perform such analyses.
Moreover, we would like to point out that linear model like the one we used has its limitation and may not be suitable for non-linear pattern discovery. We will test other non-linear models in the future. In the end, though our key driver genes have been supported by extensive literatures, some of them still need to be validated experimentally. In the future, it will be interesting to select a few less studied key driver genes for knock-out/ knock-down study in specific cell lines.

Methods
Linear regression model of gene expression changes. For each drug, we assume gene expression is affected by drug sensitivity and other confounding factors including age, sex, batch, cancer and tissue types, and the top 3 genotype principal components denoting the race and other genetic factors (see GTEx study 54 ). Similar to Yang et al. 55 we model the gene expression change using the following linear regression formula: where, Y ij is the expression level of gene j in sample i, Sensitivity i is the drug sensitivity of sample i, Age i denotes the age of sample i, Sex i denotes the sex of sample i, Batch i denotes the batch information of sample i, Cancer i denotes the cancer type of sample i, Tissue i denotes the tissue type of sample i, Genotype ik denotes the k-th principle component value of the genotype profile for the i-th sample, ε ij is the error term. In addition, β j is the regression intercept, γ j is the sensitivity regression coefficient, λ j is the age regression coefficient, μ j is the sex regression coefficient, θ j is the batch regression coefficient, κ j is the cancer type regression coefficient, φ j is the tissue regression coefficient, δ k is the regression coefficient for the k-th genotype PC. Sex i , Batch i , Cancer i , Tissue i are factor variables.
If γ j is significantly deviated from 0, gene j is considered to be drug sensitivity-associated. Gene j is positively associated with drug sensitivity if γ j > 0 and negatively associated if γ j < 0. We performed the false discovery rate (FDR) adjustment on the p-values using Benjamini Hochberg method and a FDR less than 0.1 is used as the significance threshold throughout the paper unless other thresholds are specified in special cases.

Removal of confounding factors based on principle component analysis. Removing confounding
factors is usually indispensable in revealing the true relationship between gene expression change and drug sensitivity. Besides sex, age, tissue type, batch, cancer type, major principle components (PCs) of genotype are also frequently used as confounding factors in gene expression analysis to boost true signal detection. Top genotype PCs usually covers the information of ethnicities and sub-population structures 54 . To ensure that over-fitting is not a big concern and the DRA genes we found are generally true, we performed permutation analysis. Specifically, we randomly permuted the sensitivity values of samples and repeated the test for 1,000 times. We then count the number of times that there are more DRA genes in permuted data than that in the original data and the average number of DRA genes in the permuted data, which serves as a background DRA genes detected by random.

Function enrichment analysis.
We performed functional enrichment analysis on the positively-regulated and negatively-regulated genes separately using David tools (http://david.abcc.ncifcrf.gov/summary.jsp). Benjamini score is used to control the false discovery rate and a gene set is considered to be DRA if the corresponding Benjamini score is less than or equal to 0.05.
Gender and age effect of drug sensitivity. We separated the male and female samples and analyzed DRA gene expressions using the modified regression model. Then used one tailed Fisher's exact test to calculate the significance of overlapping DRA genes between male and female on the top 100, 200, . . . , 5000 sensitivity-related genes. The sensitivity genes from male and female are considered to be significantly overlapped if the p-value of the test score is less than or equal to 0.05. Similarly, we separated the samples evenly into young and old group according to sample age and applied the modified regression model to data for separate groups respectively. As the sample size of each group is equal, we simply compared overlapping of DRA genes in each group using the Fisher's exact test.
Module Differential Connectivity (MDC) and Key Driver Analysis. We performed module differential connectivity analysis on functional modules defined by Gene Ontology 42 and Kyoto Encyclopedia of Genes and Genomes (KEGG) 43 . The genes for each GO term and KEGG pathway were retrieved using R package "org. Hs.eg.db" and "KEGG.db", respectively (on May 15, 2015). We only consider the modules with more than 30 protein coding genes. For each module, we use the R package MDC 13 to calculate the ratio of mean co-expression between the resistant and sensitive group. To evaluate the false discovery rate (FDR) for each module, we then used MDC to perform permutation tests for 50 times. Specifically, we randomly permuted the samples between the resistant and sensitive group. For each permutation, we calculated the ratio of mean co-expression between resistant and sensitive group. If MDC for original data is larger than (less than) 1, then the FDR is calculated as the number of times we have larger (smaller) MDC in the permuted group than in the original group.
The key drivers of the genes for each significant module were inferred by the package KDA 13 on protein interaction network. Specifically, the protein-protein interaction (PPI) network is constructed using the interaction confidence level 0.4 in STRING database (version 10). We first mapped the genes in a significant module into the STRING PPI network. Let A denote the set of mapped genes. For each gene in the STRING PPI network, KDA retrieved all genes with distance less than or equal to 2 from the gene. Let B denotes the gene set. A hyper geometric test is used to calculate the significance of the enrichment between gene set A and B. We also permuted genes in the 2-neighbor of the gene to calculate the FDR and the genes with FDR less than 0.05 is considered as key driver genes. We plot the key driver genes and their connectivity in the PPI network using Cytoscape (http:// www.cytoscape.org/).


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
The combined effect of PDX1, epidermal growth factor and poly-L-ornithine on human amnion epithelial cells’ differentiation

Background It has been suggested that the ectopic expression of PDX1, a dominant pancreatic transcription factor, plays a critical role in the developmental programming of the pancreas even from cells of unrelated tissues such as keratinocytes and amniotic fluid stem cells. In this study we have chosen to drive pancreatic development in human amnion epithelial cells by inducing endogenous PDX1 expression. Further, we have investigated the role of Epidermal Growth Factor (EGF) and Poly-L-Ornithine (PLO) on this differentiation process. Results Human amnion epithelial cells expressed high levels of endogenous PDX1 upon transduction with an adenoviral vector expressing murine Pdx1. Other markers of various stages of pancreatic differentiation such as NKX6.1, SOX17, RFX6, FOXA2, CFTR, NEUROD1, PAX4 and PPY were also expressed upon Pdx1 transduction. Although initial expression of pancreatic progenitor markers was higher in culture conditions lacking EGF, for a sustained and increased expression EGF was required. Culture on PLO further increased the positive impact of EGF. Conclusion Pancreatic marker expression subsequent to mPdx1 transduction suggests that this approach may facilitate the in vitro differentiation of hAECs into cells of the endocrine pancreas. This result may have important implications in diabetes therapy. Electronic supplementary material The online version of this article (doi:10.1186/s12861-016-0108-y) contains supplementary material, which is available to authorized users.

Background
The global incidence of diabetes is expected to cross the 300 million mark by 2030 [1]. With current treatment strategies being inadequate in addressing the complications associated with the disease [2], an alternative cellbased therapy approach is urgently needed. The Edmonton protocol was the first to provide a proof-of-concept of a cell-based therapy for diabetes [3]. However, due to a severe shortage of donor pancreases as well as other complications associated with the procedure itself such as the need for life-long immunosuppression and failure to achieve complete insulin independence [4], alternative approaches need to be examined. Multipotent cells from foetal tissue may provide one such alternative.
To date, attempts have been made to differentiate amnion epithelial cells (AECs), placenta-derived multipotent progenitor cells and amniotic fluid stem cells down the pancreatic lineage [5][6][7][8][9]. However, multiple protocols exist for the differentiation of the same cell type. Further, the molecular basis of the differentiation process and individual media components have not yet been studied. The present study aims to bridge this gap in current research.
The pancreatic and duodenal homeobox-1 (Pdx-1) is a master regulator of pancreatic development and differentiation [10,11]. It is responsible for the differentiation of progenitor cells into cells of the endocrine pancreas. The complex ductal architecture of the pancreas is also established by a reiterative program of branching morphogenesis which proceeds concomitantly with peak Pdx1 transcription factor expression [12]. PDX1 has been suggested to play a critical role in the developmental re-programming of the pancreas from cells of unrelated tissues such as keratinocytes and amniotic fluid stem cells [6,13]. Furthermore, induced PDX1 expression in human and non-human primate amniotic fluidderived stem cells (AFSCs) cultured on PLO-coated plates, was shown to initiate pancreatic differentiation [3,4]. PDX1 is thus, an important factor in the development of both the pancreatic islet and ductal cells.
Epidermal growth factor (EGF) is a standard component of human amnion epithelial cell growth media as it improves the growth of these cells [14,15]. In the developing pancreas, EGF increases the mass of pancreatic epithelial cells whereas a lack of it promotes endocrine differentiation [16,17]. When administered to adult pigs, excess of EGF leads to pancreatic duct hyperplasia [18]. One aspect of this study therefore is to study the effect of EGF on hAECs that are being differentiated down the pancreatic lineage.
PLO is a synthetic amino acid chain that is positively charged and widely used as a coating to enhance cell attachment and adhesion to both plastic ware and glass surfaces. PLO has previously been used in the differentiation of embryonic stem cells into insulin-producing clusters [19,20].
In this study, we have used ectopic expression of mPdx1 as a potential approach for the differentiation of hAECs into pancreatic progenitors. We found that endogenous PDX1 expression was induced several fold upon mPdx1 transduction. Several other genes that are expressed by pancreatic progenitor cells such as SOX17, RFX6, CFTR and FOXA2 were also expressed. The presence of EGF and PLO in the culture environment potentiated this expression. A two-tailed t-test was used to determine the statistical significance of the observed changes in gene dynamics.

Results
Baseline expression of pancreatic regulatory genes in hAECs (passage 2) showed high CK19 and SOX9 expression (Cq <25) (Fig. 1). In fact CK19 expression was higher in hAECs compared to adult human islets. Moderate RFX3 expression (Average Cq 27) was also observed. Expression of all other genes that were tested was either very low or absent in hAECs (Cq ≥35). We therefore attempted to initiate the process of pancreatic differentiation of hAECs by transient transduction of mPdx1, a critical factor in the development of the pancreas.

Effect of mPdx1 transduction
Transduction of hAECs with non-integrating, recombinant adenovirus harbouring the mouse Pdx1 gene was successful. Transduction efficiency increased in a doseand time-dependent manner from 10 to 200 multiplicity of infection (MOI) of the adenoviruses (Additional file 1). However, since higher concentrations of the adenoviral vector caused increased cell death over an extended culture period, we performed all subsequent experiments Fig. 1 Comparison of average Cq values between untransduced p2 hAECs and crude human islet preparations for pancreatic progenitor-specific genes. Pancreatic progenitor-specific genes were assayed by means of qPCR for untransduced p2 hAECs (baseline gene expression) and crude human islet preparations (positive control). Error bars represent mean ± SD. It can be noted that endogenous CK19 and SOX9 expression is high in untransduced hAECs while expression of PDX1, RFX6, SOX17, FOXA2 and NKX2.2 is virtually absent (Cq > 35-gene expression undetectable, Cq 30-35-low level of gene expression, Cq 25-30-high level of gene expression, Cq <25-very high level of gene expression) with 50 MOI of the virus. Transduction efficiency was 12 % at 24 h and 69 % at 48 h at this viral titre.
Subsequent to mPdx1 transduction, there was a dramatic increase in expression of endogenous human PDX1. At day 2, expression was almost 200-fold higher than untransduced controls (Fig. 2a). Expression increased marginally up to 7 days post-transduction. The expression of the PDX1 gene was confirmed by comparing PDX1 gene expression in cells that were transduced with a control EGFP adenovirus. EGFP transduction did not cause expression of PDX1 (Data not shown).
Expression of PDX1 in turn caused many other pancreatic marker genes to be expressed (Additional file 2). In particular, there was a greater than 200-fold increase in SOX17 expression and a greater than 1000-fold increase in RFX6 expression in mPdx1 transduced cells as compared to untransduced controls on day 2 posttransduction (Fig. 2a). Although there was a further 2fold increase in SOX17 expression, RFX6 expression dropped by approximately 3-fold at the end of the culture period. RFX3 expression also dropped significantly by the 7th day (Fig. 2c). Interestingly, the pancreatic endocrine progenitor markers-NGN3, NEUROD1 and PAX4 began to be expressed by day 7, albeit at low levels ( Fig. 2b). Significant expression of markers of β-cell committed cells viz. NKX6.1, NKX2.2 and PAX6 were also observed at the 2 day and 7 day time points although the difference in expression of these genes between the two time points was not significant. Further, expression of adult endocrine pancreas marker genes, GCG and PPY, also began by day 7 although this seemed to be more of an effect of the adenoviral transduction itself because even EGFP expressing cells expressed these two markers.
Microscopically, it was observed that cell proliferation ceased upon transduction with mPdx1. This was not the case with control cells transduced with EGFP (Additional file 3).
Expression of RFX6 and RFX3 at the protein level was confirmed by Western blot analysis (Additional file 4).

Effect of EGF on mPdx1 transduction
Since previous reports have suggested that an environment lacking EGF drives the pancreatic differentiation process forward, we investigated if this is indeed the case for hAECs transduced with mPdx1. Fig. 2 Effect of adenoviral transduction of human amnion epithelial cells on pancreatic marker gene expression. Fold-change expression of human pancreatic marker genes was assayed by means of qPCR (a) 2 days and (b) 7 days after mPdx1 adenoviral transduction of p2 hAECs. Fold-change is relative to untreated cells. Error bars represent mean ± SD. Only those genes with a Cq <35 are represented here. Of the 21 genes assayed, expression of 13 genes on day 2 and 18 genes on day 7 were detected. (c) Five genes in particular-NANOG, SOX17, NGN3, NEUROD1 and PAX4-showed significant increase in expression while RFX6 and RFX3 showed a significant loss of expression over the 7 day period. (* p < 0.05, ** p < 0.005, *** p < 0.0005; n = 4) In general it was observed that a lack of EGF in the culture medium, did not make any statistically significant difference to gene expression ( Fig. 3 and Additional file 5) although by day 7 there seemed to be a loss of gene expression across the marker panel, including the endocrine pancreas markers. However this was observed as a less than 1-fold change in expression in 16 of 18 genes tested. Subsequent experiments therefore continued to be carried out in the presence of 10 ng/mL EGF.

Combined effect of EGF and PLO on mPdx1 transduction
Pancreatic differentiation experiments with AFSCs have previously shown low level of insulin mRNA expression only in cells grown on a PLO coating [5]. We wanted to check therefore if growing the adenovirally transduced hAECs on PLO-coated plates would potentiate the process of pancreatic differentiation.
It was observed that the gene expression dynamics of cells growing on PLO was more or less the same irrespective of the concentration. However culture on 0.001 % PLO resulted in a statistically significant higher expression of most genes (Additional file 6). Specifically, expression of endogenous PDX1 was higher in cells cultured on 0.001 % PLO on both days as compared to their non-PLO counterparts as well as cells grown with a higher concentration of PLO (Fig. 4). Contrastingly, RFX3 expression was higher in 0.01 % PLO cultures on both days. Expression of other genes such as NKX6.1, NKX2.2, NEUROD1, PAX4, SOX17, PAX6, GCG and PPY varied with culture condition and day of culture although their expression was significantly higher than their non-PLO counterparts. In general, the combined effect of EGF and PLO caused an increase in expression of the early gene markers such as PAX4, PAX6 and SOX17. By the 7th day, the expression of these genes fell whereas those of later genes such as NGN3 and NEUROD1 increased although NGN3 expression was not statistically significant.
Under the microscope it was observed that the lower concentration of PLO was less cytotoxic (Additional file 7).

Discussion
The human amnion is a foetal membrane that surrounds the placenta and is derived from the epiblast as early as 8 days after fertilization [21]. There are 255 births worldwide every minute [22]. More than half of the placentas obtained post-partum are discarded. Since the Fig. 3 Effect of EGF on mPdx1 transduction. The expression of human pancreatic marker genes by mPdx1-transduced hAECs in a culture environment lacking EGF was assayed by means of qPCR over a 7 day period. (a) Gene expression in cells grown in the absence of EGF for 2 days (b) Gene expression in cells grown in the absence of EGF for 7 days (c) Comparison of gene expression in cells grown in the absence of EGF for 2 and 7 days. Fold change was compared to cells grown in the presence of 10 ng/mL EGF. Although it appeared that there was a loss of gene expression across the panel by the 7th day, this loss was not significant when compared to expression in cells grown in an environment containing EGF. (* p < 0.05, ** p < 0.005, *** p < 0.0005; n = 2). Error bars represent mean ± SD placenta and its surrounding membranes act as a maternal-foetal interface, it is quite likely that multipotent progenitor cells can be obtained from them. Several studies have in fact, shown the extraordinary differentiation potential of placenta derived multi progenitor cells and human amnion epithelial cells into cells of all three dermal lineages including the endocrine pancreas, neuronal cells etc. [8,9,14,23,24]. The present study has added to the knowledge obtained from such studies.
The complex development of the pancreas involves a fine play between several transcription factors, chief amongst which is the transcription factor Pdx1 [10][11][12]. In the present study, endogenous hPDX1 expression was activated. This caused cells to stop proliferating and to undergo differentiation instead.
Once PDX1 expression was activated, many other downstream genes involved in the pancreatic differentiation process were also activated, chief amongst which was the dramatic expression of RFX6. Although Rfx6 has been suggested to be an Ngn3-dependent transcription factor [25], we observed a significant, albeit low level of NGN3 expression only at the 7 day time point by which time RFX6 expression had decreased 3-fold. A significant increase in expression of other genes downstream of NGN3 such as NEUROD1 and PAX4 were also observed at this point, although the expression of these genes was very low in terms of Cq.
General endoderm marker, SOX17, and anterior endoderm marker, FOXA2, expression were increased. It was observed that as the expression of these two genes increased over the 7 day period, there was a corresponding loss of RFX6 expression. This is consistent with Pearl et. al's observation in the Xenopus system where RFX6 was expressed in two waves, the first of which occurs in the anterior endoderm, up stream of NGN3 expression [26]. However they also positioned Rfx6 up stream of Pdx1 since they observed lower gene expression of PDX1 in RFX6 knockouts. Our observation differed in that we did not see a loss of PDX1 expression as RFX6 expression reduced. In fact, PDX1 expression continued to increase over the 7 days in cultures containing EGF and 0.01 % PLO. Corroborating evidence is provided by Suzuki et al. who observed significant loss of Rfx6 expression upon Pdx1 silencing in mouse intestinal epithelial cells but did not see any significant effect of Rfx6 silencing on Pdx1 expression [27]. This suggests that some of our cells were still in the early stages of pancreatic development similar to Soyer J et al.'s observation that Rfx6 is initially expressed in the Pdx1-positive gut endoderm cells in mice and zebrafish [25].
Upon the extended culture time period, the expression of the multipotent pancreatic duct-like cell markers of FOXA2 and CFTR increased. While there was no significant expression of endocrine cell markers such as MAFA and INS, there was a greater than 10-fold increase in expression of the alpha cell marker, GCG and PP cell marker, PPY. EGF alone could sustain this expression over the entire duration of culture. Fig. 4 Combined effect of EGF and PLO on mPdx1 transduction. Two days and 7 days after mPdx1 adenoviral transduction of p2 hAECs grown in the presence of EGF and PLO, the expression of human pancreatic marker genes was assayed by means of qPCR. Expression was compared to cells that were cultured only in the presence of EGF. mRNA expressions were normalized to endogenous GAPDH expression. Similar to transduced cells cultured with EGF alone, the 7 genes-NKX2.2, PAX6, PDX1, RFX6, SOX17, CFTR and FOXA2-showed a dramatic and higher change in expression upon mPdx1 transduction. In most cases, expression of genes in cells grown on PLO was higher than their non-PLO counterpart. (* p < 0.05, ** p < 0.005, *** p < 0.0005; n = 2) Since hAECs have a high endogenous CK19 and SOX9 expression (Fig. 1), the fold-change in mRNA expression of these two genes upon adenoviral transduction even with culture on PLO, was not significant. However, we hypothesize that the high endogenous expression of CK19 and SOX9, coupled with the induced expression of PDX1 and other pancreatic marker genes, lead the cells down the pancreatic lineage. Indeed, since the cells were positive for both SOX9 and PDX1, markers of the posterior foregut, we hypothesize that some of the cells are similar to E10.5-11.5 mouse embryonic pancreatic epithelium [12] or primary pancreatic progenitor cells [28]. Interestingly, Wescott et al. observed that Pdx1 hi /insulin-negative cells may represent branching epithelia [12] and Cardinale V et al. identified multipotent stem/progenitors in high numbers in the hepato-pancreatic common duct [29]. The cells which they identified were positive for endodermal transcription factors viz. SOX9, SOX17, FOXA2, PDX1, NGN3; stem/progenitor surface markers viz. EpCAM, NCAM, CD133, CXCR4; and sometimes weakly, adult liver, bile duct, and pancreatic genes viz. albumin, CFTR and insulin, respectively. Since we were able to observe expression of some these markers in our cells as well, it provides further evidence to the possibility of our cells representing a mixed population of cells along varying degrees of pancreatic differentiation some of which are similar to the multipotent progenitors of the pancreatic duct.
Although earlier studies have indicated that a lack of EGF in the culture environment drives the differentiation process [16], there are some others which have suggested a role for EGF and its family members in establishing islet architecture [10]. In the present study we found that EGF was required to sustain and in some cases even increase the expression of pancreatic marker genes. While a lack of EGF reduced the stemness of the differentiating hAECs, it also caused a significant loss of later stage pancreatic markers including the alpha cell marker, GCG.
Human and non-human primate amniotic fluidderived stem cells (AFSCs) have previously been induced to express PDX1 when cultured on PLO-coated plates. This was shown to initiate pancreatic differentiation [5,6]. In this study, it was observed that while PLO itself did not initiate pancreatic differentiation, it greatly potentiated pancreatic differentiation upon mPdx1 transduction. What was interesting to observe was that in conditions where a lack of EGF resulted in a higher expression of a particular gene marker, a concentration of 0.01 % PLO provided similar results. When EGF or a lower concentration of PLO was introduced into the culture environment, the gene expression levels fell. It seems therefore, that a delicate balance between EGF and PLO needs to be achieved in order to drive differentiation as opposed to the maintenance of cell status.
With the incidence of diabetes on the rise and current therapy still being inadequate in preventing disease progression, the bioartificial pancreas could well be the future of diabetes therapy with enormous potential for a cure. The in vivo differentiation potential of the heterogenous population of pancreatic progenitors derived from multipotent hAECs in this present study is warranted as it could suggest their future use for cellular therapy for insulin-dependent diabetes.

Conclusion
Since the presence of true pancreatic stem cells is still a matter of debate, the plasticity of other, easily available multipotent cells could be important for a cell-therapy based approach to treat insulin-dependent diabetes. We propose human amnion epithelial cells that have been engineered to over express PDX1, as potential sources of pancreatic multipotent progenitors and thus as precursors for endocrine pancreas development. This process can be further potentiated with a mix of epidermal growth factor and poly-L-ornithine in the culture environment. This has important implications in the field of diabetes therapy.
Human pancreases were procured from deceased donors under an MOU with the Carolina Donor Services (CDS) for the supply of organs from individuals who had signed up to donate their organs for research. The CDS is the federally designated organ procurement organization, including hospitals and transplant centers that perform heart, lung, liver, kidney, pancreas and intestine transplantation. Donated organs are processed after obtaining authorization from the donor or their family.
Crude preparations of adult human islets from donor pancreases were made by a modified method of Ricordi C et al. [30]. Isolated tissue pellet was immediately frozen without any preservation media at −80°C until use for RNA extraction. The research was approved by the Institutional Review Board committee of the Wake Forest School of Medicine (#IRB00028826).

Adenoviral transduction and study on the effect of EGF
Adenovirus expressing mouse Pdx1 (mPdx1) was a gift from Drs. Christopher Newgard and Sarah Ferber at Duke University.
To confirm adenoviral transduction by fluorescence microscopy, p0 hAECs were seeded at a density of 5,000 cells/cm 2 in chamber slides (1.7 cm 2 per chamber) and cultured under standard conditions overnight. The next morning, the cultures were washed twice with DPBS or plain DMEM-LG and the appropriate viral titre diluted in 500 μL plain DMEM-LG was added. Four titres were used-10, 50, 100 and 200 MOI. After 6 h under standard culture conditions, an equal volume of complete medium was added to the cultures and cells were cultured thus overnight. The next morning, virus containing media was removed and cells were washed twice with DPBS or plain DMEM-LG. Two millilitre complete medium was added to each chamber and cells were cultured under standard conditions. The cells were fixed and immunostained 24 h and 48 h after removal of virus-containing media. Primary antibody used for mPdx1 immunofluorescence staining was a monoclonal anti-mouse Pdx1 antibody (R&D systems, USA; catalog # MAB2419). Secondary antibody was horse anti-mouse IgG-Texas red conjugated antibody (H + L; Vector Labs, USA; catalog # TI-2000). All nuclei were counter stained with DAPI (Vector Labs, USA; catalog # H-1200). Number of Texas red and DAPI stained nuclei were counted manually in two different fields. Efficiency of transduction was then estimated as: % transduction ¼ ðNumber of Texas red stained nuclei =Number of DAPI stained nucleiÞ Â 100 For RNA isolation and qPCR experiments, p2 hAECs that had achieved confluence were trypsinized, counted and plated in 6-well plates at a seeding density of 10,000 cells/cm 2 . Adenoviral transduction with 50 MOI of mPdx1-adenovirus was carried out as described above. One set of cultures lacked EGF in the complete medium. Cells were cultured for a further 2 days or 7 days in standard culture medium before RNA was extracted from them. Cells were also observed microscopically daily in order to determine if there were any accompanying morphological changes. At the time of cell culture, cultures were performed in biological triplicates.

Adenoviral transduction and study on the combined effect of EGF and PLO
One set of plates were either coated with 0.01 % or 0.001 % Poly-L-Ornithine (PLO; Sigma-Aldrich, catalog # P2533) according to the CSH protocols [31]. Cells were then cultured and transduced with mPdx1-adenovirus as described above. Cells were cultured for a further 2 days or 7 days in standard culture medium containing 10 ng/ mL EGF. Cells were also observed microscopically daily in order to determine if there were any accompanying morphological changes. At the time of cell culture, cultures were performed in biological triplicates.

RNA extraction
At the 2 day or 7 day time point, RNA was isolated from the experimental and control cells using the 5 PRIME manual PerfectPure TM RNA Cultured Cell Kit (5 PRIME, Inc., USA; catalog # 2302340). Since it was found that RNA yield from cells in one well of the 6-well plate was low, RNA from 3 wells were pooled together to represent one sample.
RNA was also isolated from p2 hAECs that had been grown in the absence of adenovirus, EGF and PLO. RNA was isolated from frozen adult human islets when required.
Total RNA content was estimated using a Nanodrop 2000c (Thermo Scientific, USA).
cDNA synthesis 500 ng of total RNA was used per 20 μL reaction, to synthesise cDNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, USA; catalog # 4368813). PCR conditions were as provided by the manufacturer's protocol i.e. 25°C for 10 min, 37°C for 2 h and 85°C for 5 min. The sequence of the primers that were obtained from Integrated DNA technologies (IDT, USA) and purified by standard desalting are given in Table 1. The specificity of the primers for the particular target gene was evaluated in silico using the Primer BLAST tool [32].

Quantitative real-time PCR (qPCR)
cDNA equivalent to 12.5 ng of initial RNA was used per 20 μL reaction for qPCR analysis. RFX6, RFX3, CFTR, CK19 and SOX9 expression were all estimated using the SYBR® Green method. 250 nanomolar of each primer (forward or reverse) was used per reaction. In the case of hydrolysis probes (  [33,34]. hGAPDH was used as the reference gene for all experiments. Unmanipulated p2 hAEC (untransduced; grown in the absence of EGF and PLO) gene expression was the quantification calibrator (baseline gene expression) for experiments in which the effect of mPdx1 transduction was studied. For all other experiments, mPdx1-transduced hAECs grown in the presence of EGF was used as the quantification calibrator. All results are represented according to the MIQE guidelines [35].
For the present study, gene expression in terms of quantification cycle (Cq) was classified as high (Cq <25), medium (Cq between 25 and 30), low (Cq between 30 and 35) and no expression (Cq ≥ 35). For each culture condition, only those genes that had a Cq <35 are represented graphically.
For individual gene expression under each culture condition, the statistical significance was estimated by a two tailed t-test (assuming equal variances) in Microsoft Excel. The 2^-Cq values for each replicate were used for the calculations. Alpha was set at 0.05.

Ethics approval and consent to participate
Human pancreases were procured from deceased donors under an MOU with the Carolina Donor Services (CDS) for the supply of organs from individuals who had signed up to donate their organs for research. The CDS is the federally designated organ procurement organization, including hospitals and transplant centers that perform heart, lung, liver, kidney, pancreas and intestine transplantation. Donated organs are processed after obtaining authorization from the donor or their family.
The use of human amnion epithelial cells was approved by the Institutional Review Board committee of

Consent for publication
Not applicable

Availability of data and materials
The data sets supporting the results of this article are included as additional files.

Additional files
Additional file 1: Figure S1. Additional file 3: Figure S2. Microscopic comparison of EGFP and mPdx1 transduced cells. Cells that were transduced with either EGFP or mPdx1 were observed 2, 4 and 7 days post-transduction. It was observed that while EGFP-expressing cells continued to proliferate over the 7 day period, proliferation ceased in cells transduced with mPdx1. Cells transduced with mPdx1 thus seemed to stop the process of proliferation and continued the process of differentiation instead. Cells were viewed using the 10X objective of an Olympus inverted fluorescence microscope. (PNG 1561 kb) Additional file 4: Figure S4. Western blot analysis. In order to prove protein expression of some of the gene expression observed by qPCR, Western blot analysis was carried out with anti-RFX6 and anti-RFX3 antibodies. (a) Comparison of the nuclear and cytoplasmic extract of adult human islets and untransduced hAECs. The islet nuclear extract (INE) and to some extent the cytoplasmic (ICE) extract stained positive for RFX6. hAEC nuclear (AECNE) and cytoplasmic (AECCE) extract were negative for RFX6 staining. This is in corroboration with qPCR analysis from Fig. 1. (b) Cell lysate from day 2 mPdx1 transduced cells grown in the presence of EGF were run on a 10 % SDS-PAGE and then immunoblotted with an anti-RFX3 or anti-RFX6 antibody. RFX3 signal was stronger than RFX6 corresponding to the higher expression of the former seen even in qPCR experiments. Anti-RFX6 antibody (Rabbit polyclonal, Sigma, catalog # SAB1402062) and anti-RFX3 antibody (Mouse polyclonal, Sigma, catalog # SAB1400241) were used at 1:1000 dilutions. Corresponding secondary antibodies used were anti-rabbit IgG-HRP conjugate (Cell signaling, catalog # 7074 s, 1:1000 dilution) and anti-mouse IgG-HRP conjugate (Cell signaling, catalog # 7076, 1:1000 dilution). The images were developed using the SuperSignal™ West Femto Maximum Sensitivity Substrate (Thermo Scientific TM , catalog # 34094). (PNG 169 kb) Additional file 5: Effect of EGF on mPdx1 transduced hAECs. The qPCR data obtained upon transducing p2 hAECs with 50MOI of the mPdx1 adenovirus is given. One set of cells were grown in the absence of EGF. The statistical calculations are also included. (XLSX 66 kb) Additional file 6: Combined effect of EGF and PLO on mPdx1 transduced hAECs. The qPCR data obtained upon transducing p2 hAECs with 50MOI of the mPdx1 adenovirus is given. One set of cells were grown in the presence of 0.01 % PLO while another was grown in the presence of 0.001 % PLO. All cells were grown in the presence of 10 ng/mL EGF. The statistical calculations are also included. (XLSX 108 kb) Additional file 7: Figure S3. Microscopic observation of cells. Cells that were transduced with mPdx1 and cultured under various culture conditions were observed at the time of transduction (0th hour) and 2 days and 7 days post-transduction. It was observed that although all cells were equally healthy at the 0th hour, the various culture conditions caused changes in their health, although not morphology, from the 2nd day post-transduction. In general, cells started senescing by the 7th day with 0.01 % PLO being the most cytotoxic of all culture conditions. Cells were viewed using the 10X objective of an Olympus inverted fluorescence microscope. (PNG 763 kb)


text:
Comparative Evaluation of Three Homogenization Methods for Isolating Middle East Respiratory Syndrome Coronavirus Nucleic Acids From Sputum Samples for Real-Time Reverse Transcription PCR.

Background Real-time reverse transcription PCR (rRT-PCR) of sputum samples is commonly used to diagnose Middle East respiratory syndrome coronavirus (MERS-CoV) infection. Owing to the difficulty of extracting RNA from sputum containing mucus, sputum homogenization is desirable prior to nucleic acid isolation. We determined optimal homogenization methods for isolating viral nucleic acids from sputum. Methods We evaluated the following three sputum-homogenization methods: proteinase K and DNase I (PK-DNase) treatment, phosphate-buffered saline (PBS) treatment, and N-acetyl-L-cysteine and sodium citrate (NALC) treatment. Sputum samples were spiked with inactivated MERS-CoV culture isolates. RNA was extracted from pretreated, spiked samples using the easyMAG system (bioMérieux, France). Extracted RNAs were then subjected to rRT-PCR for MERS-CoV diagnosis (DiaPlex Q MERS-coronavirus, SolGent, Korea). Results While analyzing 15 spiked sputum samples prepared in technical duplicate, false-negative results were obtained with five (16.7%) and four samples (13.3%), respectively, by using the PBS and NALC methods. The range of threshold cycle (Ct) values observed when detecting upE in sputum samples was 31.1–35.4 with the PK-DNase method, 34.7–39.0 with the PBS method, and 33.9–38.6 with the NALC method. Compared with the control, which were prepared by adding a one-tenth volume of 1:1,000 diluted viral culture to PBS solution, the ranges of Ct values obtained by the PBS and NALC methods differed significantly from the mean control Ct of 33.2 (both P<0.0001). Conclusions The PK-DNase method is suitable for homogenizing sputum samples prior to RNA extraction.

INTRODUCTION
irrespective of the length of time between symptom onset and sample collection [2]. Real-time reverse transcription PCR (rRT-PCR) with sputum samples is commonly used to diagnose MERS-CoV infection.
If a negative result is obtained for patients with a high index of suspicion for MERS-CoV infection, additional samples should be collected and tested [3]. A series of negative results does not rule out the possibility of MERS-CoV infection [1,3]. Possible explanations for false-negative MERS-CoV test results include low viral load very early or late in the illness, inadequate sputum samples due to poor sample collection from patients with dry cough or no cough, inappropriate sample handling or shipping, and technical issues, such as viral genome mutation or PCR inhibition [3]. Sufficient nucleic acid extraction and the removal of substances that inhibit amplification are critical factors influencing the detection of MERS-CoV by rRT-PCR [4]. However, rRT-PCR of sputum samples yields conflicting results owing to the difficulty of extracting RNA from sputum containing mucus. Many reports have indicated that lower respiratory samples such as bronchoalveolar lavage (BAL) fluid and sputum contain the highest viral loads [3,[5][6][7]. However, one report indicated that no significant difference occurred in viral load cycle threshold (Ct) values between sputa and nasopharyngeal aspirates (NPAs) [8]. For this reason, sputum sample homogenization is desirable prior to isolating nucleic acids for MERS-CoV rRT-PCR. In this study, we compared three different sputum homogenization methods and identified an optimal homogenization method prior to the nucleic acid isolation from sputum samples.

Sputum quality
Forty-eight consecutive sputum samples sent to the laboratory for MERS-CoV rRT-PCR testing were evaluated in terms of the number of white blood cells, squamous epithelial cells, and ciliated columnar epithelial cells. The sputum samples were graded by using Murray and Washington's grading system [9].

Preparation of virus-spiked sputum samples
An inactivated viral culture of MERS-CoV was provided by the Division of Respiratory Viruses (Center for Infectious Diseases, Korean National Institute of Health, Korean Center for Disease Control). The Korean National Institute of Health did not disclose the exact titer of the cultured MERS-CoV. Since other viral respiratory pathogens and bacterial pathogens should be considered during diagnosis [2], we evaluated whether the three homogenization methods have an adverse impact on the detection of other respiratory viruses. For this reason, we evaluated three homogenization methods using NPA samples that were found to contain various respiratory viruses using the Anyplex II RV16 Detection kit (Seegene, Seoul, Korea). Next, we evaluated the three homogenization methods using sputum samples spiked with a 1:10 or 1:1,000 dilution of viral culture mixed with an identical volume of sputum. Finally, we evaluated each homogenization method using virus-spiked sputum samples prepared by adding one volume of 1:100 diluted inactivated viral culture to nine volumes of clinical sample (i.e., sputum or BAL).

Sputum-homogenization methods
Sputum samples were pretreated by using one of three methods. The first method involved proteinase K and DNase I treatment (PK-DNase), as recommended by the NucliSENS easy-MAG protocol (bioMérieux SA, Marcy-l'Etoile, France). One hundred microliters of proteinase K (1 mg proteinase K/1 mL proteinase K buffer; Promega, Madison, WI, USA) was added to 100 µL of sample and incubated for 15 min at 55°C. The samples were mixed by vortexing every 5 min. Subsequently, a 20-µL of DNase I solution (1 U/µL) was added for every 100 µL of sample, followed by incubation at 37°C for 30 min. In the second method, sputum samples were mixed with an equal volume of phosphate-buffered saline (PBS). In the third method, the sputum samples were mixed with an equal volume of N-acetyl-L-cysteine (NALC) in sodium citrate solution (prepared by dissolving 0.5g NALC in each 100 mL of 1.47% sodium citrate solution). Following the second and third methods, the samples were mixed by vortexing for 1 min.

RNA extraction and MERS-CoV rRT-PCR testing
RNA was extracted as described previously [10] by using the NucliSENS easyMAG system (bioMérieux SA, Marcy-l'Etoile, France). RNA from each pretreated sample (500 µL) was eluted in a 50-μL volume. Extracted RNA was tested by rRT-PCR using Anyplex II RV16 Detection kit (Seegene) for the preliminary test, or by rRT-PCR using DiaPlexQ MERS Virus Detection kit (Sol-Gent, Daejeon, Korea) with Bio-Rad CFX96 instrument (Bio-Rad, Hercules, CA, USA). The sequences of the primers and probes used with the DiaPlexQ kit covered the same general regions of the upE and ORF1a transcripts.

Statistical analysis
Quantitative results were expressed as mean ± SD, and groups were compared by using the one-sample t-test. P values < 0.05 were considered statistically significant. MedCalc software version 16.2 (MedCalc Software bvba, Ostend, Belgium) was used for statistical analysis.
The efficacies of the three homogenization methods tested with NPA samples are summarized in Table 1. For DNA respiratory viruses (such as adenovirus and human bocavirus), PK-DNase treatment caused a slight increase in Ct values compared with PBS or NALC treatment. When we additionally tested six adenovirus-positive NPAs and five bocavirus-positive NPAs, the Ct values obtained following PBS preparation were 1.50-6.18 (median 5.42) and 1.60-7.61 (median 3.4) lower, respectively, than those obtained following PK-DNase preparation. Among the 10 samples derived from RNA respiratory viruses, the minimum Ct values resulted from PK-DNase treatment of three samples, PBS treatment of five samples, and NALC treatment of two samples.
Using identical volumes of sputum and inactivated viral culture diluted 1:10 or 1:1,000, the effects of each homogenization method on the MERS-CoV rRT-PCR results were evaluated (Table 2). Using a 1:10 dilution of inactivated viral culture, all samples treated with PK-DNase exhibited lower Ct values than those subjected to PBS or NALC treatment. The differences in Ct values between samples treated with PK-DNase or PBS were 5.4-7.8 for upE and 4.4-6.8 for ORF1a, while those for samples treated with PK-DNase or NALC were 4.7-5.8 for upE and 4.4-5.6 for ORF1a. Using the same volume of a 1:1,000 dilution of inactivated viral culture, one purulent sputum sample showed a false-negative result after PBS treatment. The Ct values of the other two sputum samples were similar following PK-DNase, PBS, and NALC treatments.
Using a one-tenth volume of a 1:1,000 dilution of inactivated viral culture (to simulate conditions comparable to those of real clinical samples), the effects of each method on the MERS-CoV rRT-PCR results were evaluated (Table 3). Among the three sputum samples, false-negative results were detected in two samples following PBS treatment and in three samples after NALC treatment (Table 3). Although the Ct values after PK-DNase treatment tended to be lower than those following PBS or  To assess statistical significance, we tested 15 additional spiked samples, which were prepared by adding a one-tenth volume of 1:1,000 diluted viral culture (Table 4). Among the 15 spiked sputum samples prepared in duplicate (30 total samples), false-negative results were detected in five samples (16.7%) using the PBS method and in four samples (13.3%) using the NALC method. The ranges of Ct values for upE from sputum samples were 31.1-35.4 using the PK-DNase method, 34.7-39.0 using the PBS method, and 33.9-38.6 using the NALC method. Compared with the control, which were prepared by adding a one-tenth volume of 1:1,000 diluted viral culture to PBS solution, the ranges of Ct values following the PBS and NALC methods significantly differed from the mean Ct value (33.2) of the control (both P < 0.0001).

DISCUSSION
For the extraction of MERS-CoV RNA from sputum sample, homogenization of sputum by PK-DNase I treatment was superior to treatment with PBS or NALC and sodium citrate. Proteinase K is a serine protease that exhibits broad cleavage specificity [11]. The enzyme cleaves peptide bonds adjacent to carboxylic group of aliphatic and aromatic amino acids and is useful for the general digestion of proteins in biological samples. More-over, proteinase K degrades any RNases present in samples, preventing RNA degradation. Viral nucleic acids are protected from DNase I degradation by the viral protein capsids, although some enveloped viruses may be lost [12]. The genomes of DNA viruses such as adenovirus and bocavirus can be affected by proteinase K and DNase I, possibly owing to the digestion of viral capsid proteins by proteinase K and subsequent degradation of genomic DNA by DNase I.
Confirmation of probable or suspected MERS-CoV infection based on clinical and epidemiological criteria can only be accomplished through laboratory testing. However, other viral respiratory pathogens and bacterial pathogens should be considered during diagnosis [3]. When using lower respiratory samples originally earmarked for MERS-CoV testing, other viral respiratory pathogens should be tested for routine diagnostic processes. Because the ability to detect DNA respiratory viruses can be reduced by DNase I treatment, proteinase K treatment alone is worth considering for the routine detection of respiratory viruses from sputum samples.
The United States Centers for Disease Control and Prevention recommends that respiratory samples be extracted by using either the NucliSENS easyMAG system or MagNA Pure Compact Nucleic Acid Isolation kit [2]. Sputum samples categorized into groups 4-5 according to Murray and Washington's grading system [9] usually cannot undergo automated extraction owing to high viscosity and, thus, require a pre-homogenization process.  Treatment with PBS or NALC and sodium citrate occasionally produces false-negative results with purulent sputum. Adding PBS and mixing via vortexing are necessary to culture respiratory viruses from sputum, and this procedure offers the advantage of involving minimum sample handling. However, thoroughly homogenizing thick or purulent sputum using this method is problematic. Because viral isolation by cell culture is not recommended for the routine confirmation of cases of MERS-CoV [2,3], the PBS plus vortexing method for processing sputum samples is not suitable for detecting MERS-CoV by rRT-PCR. Previously, it was reported that NALC and sodium citrate treatment (without NaOH) did not damage mycobacterial mRNA or rRNA in sputum samples [13]. NALC and sodium citrate solutions are readily available in clinical microbiology laboratories. Treatment with NALC and sodium citrate solution resulted in complete homogenization of sputum samples, including purulent sputum. However, this method is not suitable for detecting RNA targets in purulent sputum. Because proteinase K treatment requires heat inactivation of potential pathogens (incubation for 15 min at 55°C), this method offers additional laboratory safety. Only 10% of sputum samples sent for MERS-CoV rRT-PCR testing were clinically relevant in terms of bacterial culture (group 5 according to Murray and Washington's grading system [9]). Among patients infected with MERS-CoV, a dry cough is more common than a productive cough [16,17]. Therefore, optimum sample collection is not feasible for patients unable to produce sputum. However, two thirds of sputum samples in our study contained identifiable ciliated columnar epithelial cells. Dipeptidyl peptidase 4, the human receptor for MERS-CoV [16], is expressed in type I and II alveolar cells, ciliated or non-ciliated bronchial epithelial cells, and bronchial submucosal glands [18,19]. Although it is difficult to obtain high-quality sputum from patients with no cough or a dry cough, obtaining and testing sputum samples containing ciliated columnar epithelial cells is recommended, if possible.
Our study has several limitations. First, it was performed by using spiked sputum samples; thus, our results may differ from those obtained by using authentic MERS-CoV-positive samples. However, the spiked sputum samples were similar to real sputum samples, as they were produced by adding a one-tenth volume of spiked viral solution to total sputum samples. Second, we tested a relatively small number of sputum samples under each test condition. Therefore, our findings should be confirmed by additional, larger studies.
In conclusion, the PK-DNase I treatment is suitable for the homogenization of sputum prior to RNA extraction for MERS-CoV rRT-PCR detection. Furthermore, adequate sputum homogenization prior to the isolation of nucleic acids is critical for accurate diagnosis of MERS-CoV infection.