Prediction of regulatory targets of alternative isoforms of the epidermal growth factor receptor in a glioblastoma cell line

Background The epidermal growth factor receptor (EGFR) is a major regulator of proliferation in tumor cells. Elevated expression levels of EGFR are associated with prognosis and clinical outcomes of patients in a variety of tumor types. There are at least four splice variants of the mRNA encoding four protein isoforms of EGFR in humans, named I through IV. EGFR isoform I is the full-length protein, whereas isoforms II-IV are shorter protein isoforms. Nevertheless, all EGFR isoforms bind the epidermal growth factor (EGF). Although EGFR is an essential target of long-established and successful tumor therapeutics, the exact function and biomarker potential of alternative EGFR isoforms II-IV are unclear, motivating more in-depth analyses. Hence, we analyzed transcriptome data from glioblastoma cell line SF767 to predict target genes regulated by EGFR isoforms II-IV, but not by EGFR isoform I nor other receptors such as HER2, HER3, or HER4. Results We analyzed the differential expression of potential target genes in a glioblastoma cell line in two nested RNAi experimental conditions and one negative control, contrasting expression with EGF stimulation against expression without EGF stimulation. In one RNAi experiment, we selectively knocked down EGFR splice variant I, while in the other we knocked down all four EGFR splice variants, so the associated effects of EGFR II-IV knock-down can only be inferred indirectly. For this type of nested experimental design, we developed a two-step bioinformatics approach based on the Bayesian Information Criterion for predicting putative target genes of EGFR isoforms II-IV. Finally, we experimentally validated a set of six putative target genes, and we found that qPCR validations confirmed the predictions in all cases. Conclusions By performing RNAi experiments for three poorly investigated EGFR isoforms, we were able to successfully predict 1140 putative target genes specifically regulated by EGFR isoforms II-IV using the developed Bayesian Gene Selection Criterion (BGSC) approach. This approach is easily utilizable for the analysis of data of other nested experimental designs, and we provide an implementation in R that is easily adaptable to similar data or experimental designs together with all raw datasets used in this study in the BGSC repository, https://github.com/GrosseLab/BGSC. Electronic supplementary material The online version of this article (10.1186/s12859-019-2944-9) contains supplementary material, which is available to authorized users.


Background
Glioblastoma is the most malignant and most frequent primary cerebral tumor in adults and is responsible for 65% of all brain tumors [1]. One potential molecular target amplified in 36% of glioblastoma patients is the epidermal growth factor receptor (EGFR), and the expression of EGFR is associated with prognosis in cancer [2]. EGFR is known to affect growth and survival signals and to play a crucial role in the regulation of cell proliferation, differentiation, and migration of various tumor entities [3]. Hence, EGFR is well known as a prognostic tumor marker and therapeutic target in different tumor entities.
The full-length transmembrane glycoprotein isoform of EGFR consists of three functional domains of which the extracellular domain is capable of binding at least seven different ligands such as EGF, AREG, or TGF-α [4]. However, there are at least three different truncated EGFR splice variants (II, III, and IV). Up to now, only the fulllength EGFR isoform I translated from EGFR splice variant I is well investigated, but comparatively little is known about the biological significance of the truncated EGFR isoforms II-IV translated from EGFR splice variants II-IV.
EGFR isoforms II-IV lack the intra-cellular tyrosinekinase domain [5], and Maramotti et al. [6] describes that EGFR isoforms II-IV can potentially function as natural inhibitors of EGFR isoform I. EGFR isoforms II-IV bind EGF with similar binding kinetics but lower binding affinity than EGFR isoform I [7], which binds EGF with a dissociation constant of 1.77 × 10 −7 M [8].
Different tumor therapies targeting EGFR via antibodies or small molecules often do not have response rates as successful as expected. EGFR isoforms II-IV may be responsible for therapeutic failures because they do not contain the tyrosine-kinase domain targeted by small molecules. However, they do contain the extracellular N-terminus of EGFR, which is bound by therapeutic antibodies. Nevertheless, EGFR-specific antibody therapy requires the interaction of EGFR-bound therapeutic antibodies with presenting cells. EGFR isoforms II-IV are soluble proteins that do not mark the expressing cell itself, but rather diffuse in the extracellular space, probably bind to surrounding non-tumor cells, and possibly mislead the immune system. This problem motivated the present work of perturbing the profile of the four EGFR splice variants using small interfering RNAs (siRNAs) that differentially target these splice variants and of measuring the resulting expression responses using traditional microarrays. It is impossible to knock-down only EGFR splice variants II-IV and not EGFR splice variant I by RNAi because there is no region specific to only EGFR splice variants II-IV. Hence, we performed the RNAi experiments according to the nested experimental design as shown in Table 1. RNAi by siRNA ALL x 5 x 6 The six corresponding logarithmic expression values per gene are denoted by Based on this design, the associated effects of a knockdown of EGFR splice variants II-IV can only be inferred indirectly by subtracting the effects found by knocking down only EGFR splice variant I from the effects found by knocking down all EGFR splice variants I-IV. The problem of only indirectly measurable gene regulation or receptor effects of nested splice variants is widespread in many regulatory pathways and many species, so we developed a two-step bioinformatics approach for the prediction of putative target genes called Bayesian Gene Selection Criterion (BGSC) approach, which we tested by quantitative real-time polymerase chain reaction (qPCR) experiments. The rest of this paper is structured as follows: In Results, we describe the identification of a cell line with an inducible EGFR-signaling pathway, investigate the specificity of siRNAs, introduce the two-step BGSC approach for predicting putative target genes regulated by EGF via EGFR isoforms II-IV and not by the full-length EGFR isoform I or other receptors, and describe the qPCR validation experiments. In Discussion, we discuss the adjustability of the EGFR-signaling pathway in cell line SF767 and the biological relevance of the validated genes.

Identification of a cell line with an inducible EGFR-signaling pathway
A meaningful analysis of the EGFR-signaling pathway is possible only in a cell line with an adjustable pathway, e.g., by a response to ligand stimulation or treatment by a tyrosine kinase inhibitor (TKI) [9]. Hence, we investigated four glioblastoma cell lines in a pilot study to identify a cell line with an adjustable EGFR-signaling pathway. Figure 1 shows the measured protein levels of phosphorylated AKT (pAKT) resulting from the treatment of two of these cell lines U251MG and SF767 with increasing levels of recombinant ligand EGF. We found that the pAKT (Ser473) level in cell line U251MG is constantly high, possibly resulting from the mutated PTEN gene [10]. In the PTEN wild-type cell line SF767 [11], pAKT showed a level of activity even without adding recombinant EGF due to the E545K-mutation of gene PIK3CA present in this cell line [12]. However, the activity of pAKT could be Fig. 1 Western blot analysis of the two glioblastoma cell lines U251MG and SF767. U251MG is a PTEN mutant and PIK3CA wild-type cell line and SF767 is a PTEN wild-type and PIK3CA (E545K) mutant cell line. Cells were treated for 24 hours with different levels of the EGFR-ligand EGF (0-50 ng/ml). The levels of HER2 and EGFR are reduced by EGF-dependent degradation of the formed and internalized EGF-HER2/EGFR complexes. The activation of AKT-protein (phosphorylation of the Ser473) is detectable in an EGF-dependent manner in cell line SF767, whereas the pAKT level is constantly high in cell line U251MG. These observations indicate that the EGFR-signaling pathway is inducible in cell line SF767, but not in cell line U251MG. Anti-β-actin staining was done as a loading control, and BIRC5 (survivin) was used as an indicator for proliferation activity increased three-fold by adding recombinant EGF as a ligand, indicating that the EGFR-AKT signaling pathway was inducible in an EGF-dependent manner (Fig. 1). Figure 1 also shows that the full-length EGFR protein disappeared by applying a high concentration of EGF of 50 ng/ml to cell line SF767. This high concentration of EGF leads to the saturation of the full-length EGFR protein with the ligand EGF, to the subsequent internalization and degradation of the formed EGF-EGFR complex, and thus to the observed disappearance of the full-length EGFR protein.

Specificity of siRNAs
We performed RNAi experiments with a siRNA against EGFR splice variant I, henceforth called siRNA I and with a siRNA against all EGFR splice variants, henceforth called siRNA ALL (Table 2). To investigate the specificity of the two siRNA constructs siRNA ALL and siRNA I , we analyzed mRNA levels and protein levels of EGFR. Figure 2 shows that the treatment of SF767 cells with the two siRNAs reduced the level of full-length EGFR protein 24 hours and 48 hours after the start of the experiment. We then analyzed the siRNA-specificity by qPCR experiments for (a) all EGFR splice variants together, (b) EGFR splice variant I (full-length), (c) EGFR splice variant IV, and (d) the two genes MMP2 and GAPDH as a control. Additional file 1: Figure S.1 shows that the application of siRNA ALL and siRNA I reduced the levels of all EGFR splice variants by 70.9% on average and the levels of the full-length EGFR splice variant I by 78.1% on average. Additional file 1: Figure S.1 also shows that the application of siRNA ALL reduced the levels of EGFR splice variant IV by 69.9% on average, that the application of siRNA I did not reduce the levels of EGFR splice variant IV, and that the application of siRNA ALL and siRNA I did not reduce the levels of the two control genes.

First step of the BGSC approach -grouping of genes
The binding affinities of the three EGFR isoforms II-IV to EGF are lower than that of the full-length EGFR isoform I [7] and probably different from each other, but yet very high [7], so we assume that the high concentration of EGF of 50 ng/ml leads to the saturation of all EGFR isoforms irrespective of their different binding affinities to EGF. Hence, we make the simplifying Next, we consider for each RNAi treatment if the genes of each group would be differentially regulated after EGFstimulation. To conceptually analyze the gene expression of each group we denote by "1" a theoretical regulation (up or down) of the group after addition of EGF and denote by "0" no regulation. Further, we define groups as regulated after EGF-stimulation if there is at least one incoming edge to the group in the graphical representation (Fig. 4), and we define groups with no incoming edge as unregulated. We consider three experimental manipulations with RNAi: negative control without RNA interference, RNAi with siRNA against EGFR splice variant I, henceforth called siRNA I , and RNAi with siRNA against all EGFR splice variants, henceforth called siRNA ALL (Fig. 4).
First, we consider the negative control without RNA interference (Fig. 4a). Here, none of the EGFR splice variants are down-regulated by a siRNA, so all target genes of EGFR isoforms and target genes of other EGF receptors can be induced by EGF. Hence, we expect differential expression under EGF stimulation of genes belonging to groups B -H on the one hand and no differential expression of genes belonging to group A on the other hand.
Second, we consider RNAi treatment with siRNA I (Fig. 4b). Here, only EGFR splice variant I is downregulated by siRNA I , so only target genes of EGFR isoforms II-IV and target genes of other EGF receptors can be induced by EGF. Hence, we expect differential expression by EGF treatment of genes belonging to groups B, C, and E -H on the one hand and no differential expression of genes belonging to groups A and D on the other hand.
Third, we consider RNAi treatment with siRNA ALL (Fig. 4c). Here, all four EGFR splice variants are downregulated by siRNA ALL , so only target genes of other EGF receptors can be induced by EGF. Hence, we expect differential expression by EGF treatment of genes belonging to groups B and E -G on the one hand and no differential expression of genes belonging to groups A, C, D, and H on the other hand. Figure 5 summarizes the different expression patterns of Fig. 4. We find that the eight gene groups show only four different expression patterns, so we reduce the eight gene groups A -H to the four simplified gene groups ad, where group A becomes group a, the union of the groups Fig. 5 Reduction of the conceptual gene groups. Genes of group A are never differentially expressed by EGF treatment. Genes of group B and E -G are always differentially expressed by EGF treatment. Genes of group C and H are differentially expressed by EGF treatment in case of control treatment (no RNAi) or simultaneous treatment with siRNA I , whereas not differentially expressed by EGF treatment in case of simultaneous treatment with siRNA ALL . Genes of group D are differentially expressed by EGF treatment in case of control treatment (no RNAi), whereas not differentially expressed by EGF treatment in case of simultaneous treatment with siRNA I or siRNA ALL . We find that the eight gene groups show only four different expression patterns, so we reduce the eight gene groups A -H to the four simplified gene groups a -d, where group A becomes group a, the union of the groups B and E -G becomes group b, the union of the groups C and H becomes group c, and group D becomes group d B and E -G becomes group b, the union of the groups C and H becomes group c, and group D becomes group d.
These simplified gene groups can be easily interpreted as follows: Genes of group a are not regulated by EGF, whereas genes of groups b − d are regulated by EGF. Genes of group b are regulated by EGF only through other receptors besides EGFR isoforms. Genes of group c are regulated by EGFR isoforms II-IV and not by other receptors. And genes of group d are regulated by EGFR isoform I and not by EGFR isoforms II-IV or other receptors. Based on this reduction, we can now formulate the goal of this work as the prediction of putative target genes regulated by EGFR isoforms II-IV and not by other receptors or, more crisply, as the goal of predicting genes of group c.

Second step of the BGSC approach -classification of genes
In the second step, we classify each potential target gene into one of the four simplified gene groups z ∈ {a, b, c, d} based on the Bayesian Information Criterion, and thereby predict target genes regulated by EGF via EGFR isoforms II-IV as those classified into group c.
In this step, we apply the oversimplified, but commonly accepted, assumption that the log-transformed expression of each gene is normally distributed [13] with a gene-specific and treatment-specific mean and variance.
For each gene, we additionally assume heteroscedasticity, i.e., equality of the six variances, of the six normally distributed logarithmic expression values under each of the six experimental conditions, an assumption commonly made in the t-test, the analysis of variance, or other statistical tests. We further assume that the six means of these six normal distributions are group specific as shown in Fig. 6.
First, we assume genes of group a (not regulated by EGF) to show no differential expression under each of the b a d c  (Table 1), as manifested by equality of the six means of the six normal distributions (Fig. 5, yellow column).
Second, we assume genes of group b (regulated by EGF through other receptors besides any EGFR isoform) to show differential expression under EGF-stimulation, irrespective of RNAi treatment targeting any EGFR isoform (Fig. 5, blue column). Hence, we assume genes of group b to have two different mean logarithmic expression levels, one in samples 1, 3, and 5, and another potentially different one in samples 2, 4, and 6 ( Table 1). We denote these two mean logarithmic expression levels by μ b0 (Fig. 6b red) and μ b1 (Fig. 6b blue) respectively.
Third, we assume genes of group c (regulated by EGFR isoform II-IV and not by other receptors) to show differential expression between the negative control and siRNA ALL treatments (Fig. 5, red column) under EGFstimulation. Hence, we assume genes of group c to have two different mean logarithmic expression levels, one in samples 1, 3, 5, and 6, and another potentially different one in samples 2 and 4 (Table 1). We denote these two mean logarithmic expression levels by μ c0 (Fig. 6c red) and μ c1 (Fig. 6c blue) respectively.
Fourth, we assume genes of group d (regulated by EGFR isoform I only) to show differential expression between the negative control and siRNA I treatment (Fig. 5, green column) under EGF-stimulation. Hence, we assume genes of group d to have two different mean logarithmic expression levels, one in samples 1, 3, 4, 5, and 6, and another potentially different one in sample 2 (Table 1). We denote these two mean logarithmic expression levels by μ d0 (Fig. 6d red) and μ d1 (Fig. 6d blue) respectively.
For genes of group a we denote the two model parameters μ a and σ a of the six normal distributions by θ a = (μ a , σ a ), and for each of the three groupsz ∈ {b, c, d} we denote the three model parameters μz 0 , μz 1 , and σz of the six normal distributions by θz = (μz 0 , μz 1 , σz).
Assuming conditional independence of the six logarithmic expression levels given group z and model parameters θ z , we can write the likelihood p(x|z, θ z ) of data x given group z and model parameters θ z as a product of six univariate normal distributions with the corresponding mean μ a , or means μz 0 and μz 1 , and the corresponding variance σ 2 z (Eqs. 1 and 2). Using the maximum likelihood principle, we obtain the estimates of model parameters θ a by Eqs. 8a and 8b and of model parameters θz forz ∈ {b, c, d} by Eqs. 8c, 8d and 8e.
To illustrate this approach, we show the six measured logarithmic expression levels together with the univariate normal probability density estimated for group a and the three pairs of univariate normal probability densities estimated for each of the three groupsz ∈ {b, c, d} for gene TPR in Fig. 7. Visually, it is easy to see that the model of group c fits best the expression profile of this gene, as it yields the best separation between the two estimated means and the smallest estimated pooled variance. Consistent with this visual observation, the four corresponding likelihoods of the six measured logarithmic expression levels are p(x|a, θ a ) = 0.004, p(x|b, θ b ) = 0.035, p(x|c, θ c ) = 4.22, and p(x|d, θ d ) = 0.012, i.e., the likelihood of the six measured logarithmic expression levels of gene TPR is highest for group c.
However, performing classification through model selection based on maximizing the likelihood is problematic when the number of free model parameters is not identical among all models under comparison. In the BGSC approach, model a has two free model parameters, while models b, c, and d have three free model parameters. Hence, a simple classification based on maximizing the likelihood would give a spurious advantage to models b, c, and d with three free model parameters over model a with only two free model parameters. To eliminate that spurious advantage, we compute marginal likelihoods p(x|z) using the approximation of Schwarz et al. [14] commonly referred to as Bayesian Information Criterion (section "Probabilistic modeling of gene expression"). Applying this approximation to gene TPR we obtain the four marginal likelihoods of the six measured logarithmic expression levels p(x|a) = 0.001, p(x|b) = 0.002, p(x|c) = 0.287, and p(x|d) = 0.001. We find that the marginal likelihood for group c is highest, which is consistent with the visual observation of Fig. 7.
To obtain the approximate posterior probability p(z|x), we now simply use Bayes' formula p(z|x) = (p(x|z)p(z))/p(x) for group z ∈ {a, b, c, d}, where p(z) is the prior probability of group z, and the denominator p(x) is the sum of the four numerators p(x|z)p(z) for z ∈ {a, b, c, d}. We assume that 70% of all genes are not regulated by EGF, so we define the prior probability for group a by p(a) = 0.70, and we further assume that the remaining 30% of the genes fall equally in groups with EGF-regulation, so we define the prior probabilities for groups b, c, and d by p(b) = p(c) = p(d) = 0.1. Using these prior probabilities, we obtain for gene TPR the four approximate posterior probabilities p(a|x) = 0.016, p(b|x) = 0.008, p(c|x) = 0.973, and p(d|x) = 0.003. We find that the approximate posterior probability for group c is highest, so we finally assign gene TPR to group c.
By applying this approach of computing the four approximate posterior probabilities for each gene and assigning each gene to that group z with the highest approximate posterior probability, we classify 8449 genes to group a, 3822 genes to group b, 3143 genes to group c, and 1328 genes to group d.

Prediction of genes belonging to simplified gene group c
For simplified gene group c, we define the subset of the 1140 genes with an approximate posterior probability a b c d Fig. 7 Probability density plot of the normal distributions of TRP. For group a we mark the logarithmic expression values x 1 , . . . , x 6 of TPR with black points, which are colored according to Fig. 6a, and assume that all six logarithmic expression levels stem from the same normal distribution. In black, we plot the probability density of this normal distribution with mean and standard deviation equal to μ and σ of the six logarithmic expression levels. For groups b -d we assume that all six logarithmic expression levels stem from a mixture of two normal distributions with independent means μ 0 and μ 1 and one pooled standard deviation σ . We mark the logarithmic expression values x 1 , . . . , x 6 of TPR with points which are colored according to indicator variables from Fig. 6 g = 0 in red and g = 1 in blue and we plot the probability densities of the two normal distributions in red and blue, respectively. For group b we assume that the logarithmic expression levels x 1 , x 3 , and x 5 stem from the normal distribution with mean μ 0 (red) and x 2 , x 4 , and x 6 from the normal distribution with mean μ 1 (blue). For gene group c we assume that the logarithmic expression levels x 1 , x 3 , x 5 , and x 6 stem from the normal distribution with mean μ 0 (red) and x 2 and x 4 from the normal distribution with mean μ 1 (blue). For group d we assume that the logarithmic expression levels x 1 , x 3 , x 4 , x 5 , and x 6 stem from the normal distribution with mean μ 0 (red) and x 2 stem from the normal distribution with mean μ 1 (blue) p(c|x) exceeding 0.75 as putative target genes regulated by EGFR isoforms II-IV and not by other receptors (Additional file 2: Table S.1), and we scrutinize six of these genes in the following section. Three of these genes (CKAP2L, ROCK1, and TPR) are up-regulated with a log 2fold changeμ c1 −μ c0 > 0.5 and three of these genes (ALDH4A1, CLCA2, and GALNS) are down-regulated with a log 2 -fold changeμ c1 −μ c0 < −0.5.
To validate the 36 logarithmic expression levels x 1 , . . . , x 6 of the six genes CKAP2L, ROCK1, TPR, ALDH4A1, CLCA2, and GALNS, we perform 108 qPCR experiments comprising three biological replicates for each gene and each treatment. Figure 8 shows the 12 log 2fold changesμ c1 −μ c0 of the microarray experiments and of the qPCR experiments. We find that the six log 2 -fold changes of the microarray experiments and those of the qPCR experiments are not identical, but in good agreement, yielding a Pearson correlation coefficient of 0.99. Moreover, the error bars, computed by using the Satterthwaite approximation, of all six genes overlap between microarray experiments and qPCR experiments.
To investigate the degree to which the expression levels of these genes respond to EGF in another glioblastoma cell line, we perform triplicated qPCR experiments in the glioblastoma cell line LNZ308 with and without EGF treatment. As CLCA2 is not sufficiently expressed in cell line LNZ308 with a log-expression of −5.8 in the Cancer Cell Line Encyclopedia data [10], we stimulate cell lines SF767 and LNZ308 with EGF (50 ng/ml for 24 hours) and measure the expression of the five remaining genes by qPCR experiments. We find that the log 2 -fold changes are not identical, but in good agreement, between the two cell lines for the four genes CKAP2L, ROCK1, TPR, and GALNS, whereas they are different between the two cell lines for gene ALDH4A1 (Additional file 1: Figure S.2).

Adjustability of the EGFR-signaling pathway in cell line SF767
To analyze the function of the soluble EGFR (sEGFR) isoforms II-IV it is essential to use a cell line with an adjustable EGFR-signaling pathway. As shown in Fig. 1, the EGFR-signaling pathway is adjustable in cell line SF767 with respect to recombinant EGF stimulation, even though cell line SF767 has a PIK3CA (E545K) mutation resulting in a baseline level of AKT activation [15]. This mutation occurs in about 30% of human breast cancers, where it leads to gain-of-function mutations in gene PIK3CA that activate the PI3K-AKT-signaling pathway constantly, thereby uncoupling the EGFR response from AKT signaling [16]. However, in cell line SF767 the level of pAKT can be increased nearly three-fold in an EGF-dependent manner (Fig. 1) consistent with the observation of Sun et al. [17].
It has been suggested that glioblastoma cell lines with helical domain mutations are still sensitive to dual PI3Ki/MEKi treatment [9], which is consistent with our observation that the EGFR-signaling pathway is adjustable in cell line SF767. Also, it has been found that Gefitinib inhibited EGFR phosphorylation in U251MG and SF767 cells, whereas Gefitinib inhibited AKT phosphorylation only in SF767 cells but not in U251MG cells [18], consistent to Fig. 1. Other EGF-induced signaling pathways such as the PLCγ -signaling pathway appear to be intact in cell line SF767 too [19].
Next, we perform western blot experiments and find that both siRNAs reduce the levels of the full-length EGFR proteins (Fig. 2). By qPCR experiments we find that siRNA ALL is capable of knocking down all EGFR splice variants and that siRNA I is capable of selectively knocking down EGFR splice variant I (Additional file 1: Figure S.1). More precisely we detect a reduction by 70.9% on average for all EGFR splice variants and a reduction by 78.1% on average for EGFR splice variant I for siRNA ALL as well as for siRNA I (Additional file 1: Figure S.1). Based on similar reductions, it appears that EGFR splice variant I is the dominant splice variant. As expected, the level of EGFR splice variant IV was reduced only by siRNA ALL .

Biological context of genes predicted to belong to simplified gene group c
Next, we investigate the biological context of the six genes predicted to belong to simplified gene group c by applying the BGSC approach under the simplifying assumption of neglecting the different binding affinities of the EGFR isoforms to EGF.
The 'Cytoskeleton Associated Protein 2 Like' (CKAP2L) protein is localized on microtubules of the spindle pole throughout metaphase to telophase in wild-type cells [20], and a knock-down of CKAP2L has been found to suppresses migration, invasion, and proliferation in lung adenocarcinoma [21].
The 'Rho-Associated Protein Kinase 1' (ROCK1) is known to play an important role in the EGF-induced formation of stress fibers in keratinocyte [22] and to be involved in the cofilin pathway in breast cancer [23]. Besides, ROCK1 has been found to promote migration, metastasis, and invasion of tumor cells and also to facilitate morphological cell shape transformations through modifications of the actinomyosin cytoskeleton [24].
Depletion of the mRNA of the 'Tumor Potentiating Region' (TPR) gene by RNAi triggers G0-G1 arrest, and TPR depletion plays a role in controlling cellular senescence [25]. Also, TPR regulates the nuclear export of unspliced RNA and participates in processing and degradation of aberrant mRNAs [26], a mechanism considered important for the regulation of genes and their deregulation in cancer cells.
The ' Aldehyde Dehydrogenase 4 Family Member A1' (ALDH4A1) gene contains a potential p53 binding sequence in intron 1, and p53 is often mutated in tumor cells [27]. Moreover, ALDH4A1 was induced in a tumor cell line in response to DNA damage in a p53-dependent manner [27], and depletion of the mRNA of ALDH4A1 by siRNA results in severe inhibition of cell growth in HepG2 cells [28].
A second gene that is transcriptionally regulated by DNA damage in a p53-dependent manner is the 'Chloride Channel Accessory 2' (CLCA2) gene. Inhibition of CLCA2 stimulates cancer cell migration and invasion [29]. Furthermore, CLCA2 could be a marker of epithelial differentiation, and knock-down of CLCA2 causes cell overgrowth as well as enhanced migration and invasion. These changes are accompanied by down-regulation of E-cadherin and up-regulation of vimentin, and loss of CLCA2 may promote metastasis [29]. Also, loss of breast epithelial marker CLCA2 has been reported to promote an epithelial-to-mesenchymal transition and to indicate a higher risk of metastasis [30].
For the 'Galactosamine (N-Acetyl)-6-Sulfatase' (GALNS) gene an effect of 17β-estradiol on the expression of GALNS could be detected by qPCR experiments in a breast cancer cell line, which is a hint to a tumor association of GALNS [31].
Up-regulation of ROCK1 and TPR and down-regulation of ALDH4A1 and CLCA2 (Fig. 8) are positively associated with the processes of migration, metastasis, and invasion of tumor cells and negatively associated with proliferation. The up-regulation of CKAP2L [32] by EGFR II-IV isoforms indicates a potential link to processes of cell-cycle progression of stem cells or progenitor cells. Overall, our interpretation of the impact of EGFR isoforms II-IV on four of six validated gene transcripts is that it seems likely that these isoforms are involved in processes of migration and metastasis of clonogenic (stem) cells, which is strongly associated with a more aggressive tumor and a worse prognosis of tumor disease.
We found that the BGSC approach was capable of detecting genes putatively regulated by EGFR isoforms II-IV and not by other receptors such as HER2, HER3, or HER4 [33], so we find it tempting to conjecture that the BGSC approach could be useful for the analysis of similarly-structured data of other nested experimental designs.

Conclusions
We have performed RNAi experiments to analyze the expression of three poorly investigated isoforms II-IV of the epidermal growth factor receptor in glioblastoma cell line SF767 with an adjustable EGFR-signaling pathway, and we have developed the Bayesian Gene Selection Criterion (BGSC) approach for the prediction of putative target genes of these EGFR isoforms under the simplifying assumption of neglecting the different binding affinities of the EGFR isoforms to EGF. We have predicted 3143 putative target genes, out of which 1140 genes have an approximate posterior probability greater than 0.75, and we have tested six of these genes by triplicated qPCR experiments. These six genes include ROCK1, which is known to be associated with EGFR regulation, as well as CKAP2L, TPR, ALDH4A1, CLCA2, and GALNS. We have found that the six log 2 -fold changes of the microarray expression levels and those of the qPCR expression levels are highly correlated with a Pearson correlation coefficient of 0.99 (p-value = 0.00002), suggesting that the set of 1140 genes might contain some further putative target genes of EGFR isoforms II-IV in tumor cells. As suggested by our anonymous reviewers we like to point out that, in addition to RNAi, CRISPR/Cas knockout [34] and replacement with each isoform would be a promising strategy to discover additional functions of the soluble EGFR isoforms besides the ones described by Maramotti et al. [6]. The analysis of isoform-specific effects in combination with RNAi treatments are an elegant way to directly downregulate specific mRNA splice variants, but that often leads to a nested experimental design for which generally no standard procedure exists. The two-step BGSC procedure of first defining easily interpretable conceptual groups of genes associated with different EGFR isoforms and subsequently classifying genes based on the approximated posterior probability to these groups seems to be a promising approach in such a situation, and this approach is readily adaptable to other and more complex experimental designs. The datasets analyzed during the current study and the R-scripts for reproducing the results and plots of this work are available in the BGSC repository, https://github.com/GrosseLab/BGSC.

Western blot and qPCR analyses
Cells were treated in lysis buffer, the protein concentration was determined using the Bradford method, and western blot analysis was performed as described in [35]. Antibodies directed against EGFR (Clone D38B1), HER2/ErbB2 (29D8), and phosphoserine 473 AKT (clone D9E) were obtained from Cell Signaling Technology Inc. (Signaling, Danvers, MA, USA), antibodies directed against β-actin were obtained from Sigma (Steinheim, Germany), and BIRC5 (Survivin) antibodies (clone AF886) were obtained from R&D systems (Richmond, CA, USA). qPCR experiments were performed as described in [35]. The primer sequences are listed in Table 3.  [37] was used to read the 47317 microarray probes into R. The neqc function of limma was used to perform a background correction followed by quantile normalization, using negative control probes for background correction and both negative and positive controls for normalization. The 16,742 array probes corresponding to 14,389 genes, which displayed a significant hybridization signal (Illumina signal detection statistic at P < 0.05) in all probes were used for further analysis.

Experimental design
For investigating which genes are activated by the four EGFR isoforms I -IV in glioblastoma cell line SF767 we use RNAi, as described in section "RNAi", for a selective down-regulation of EGFR splice variants (Table 1 rows) with and without EGF treatment (Table 1 columns). Specifically, we applied the three different RNAi treatments -(i) control without RNAi, (ii) RNAi with siRNA I , and (iii) RNAi with siRNA ALL -to glioblastoma cell line SF767.
In case (i) we performed a control experiment without RNAi treatment (Table 1, first row). Here, EGFR is not down-regulated by an siRNA, so target genes of all EGFR splice variants and other EGF receptors should be differentially expressed in columns 1 and 2, i.e., they should have different logarithmic expression levels x 1 and x 2 .
In case (ii) we performed an RNAi with siRNA I , which can bind only to the full-length EGFR splice variant I ( Table 1, second row). Hence, siRNA I down-regulates splice variant I, but not the other splice variants II-IV, and in this case target genes of EGFR isoforms II-IV and of other EGF receptors should be differentially expressed in columns 1 and 2, i.e., they should have different logarithmic expression levels x 3 and x 4 .
In case (iii) we performed an RNAi with siRNA ALL , which can bind to all four EGFR splice variants, and subsequently down-regulates all four splice variants (Table 1, third row). Here, only target genes of other EGF receptors should be differentially expressed in columns 1 and 2, i.e., they should have different logarithmic expression levels x 5 and x 6 .

Probabilistic modeling of gene expression
We propose a probabilistic model for the logarithmic expression pattern x = (x 1 , . . . , x 6 ) for each of the four groups z ∈ {a, b, c, d} defined in section "First step of the BGSC approach -grouping of genes".
First, we assume that the three logarithmic expression levels x 1 , x 3 , and x 5 corresponding to no EGF treatment are similar to each other, which corresponds to the assumption that the RNAi treatment should have no effect in case of no EGF treatment. Second, we assume that the three logarithmic expression levels x 2 , x 4 , and x 6 follow the expression patterns described in section "First step of the BGSC approach -grouping of genes" and summarized in Fig. 5.
In order to mathematically formulate the model assumptions, we introduce six indicator variables g 1 , . . . , g 6 for the groupsz ∈ {b, c, d} that indicate if the six logarithmic expression levels x 1 , . . . , x 6 are expected to be different from x 1 . Specifically, we define g n = 1 if x n is expected to be different from x 1 for n = 1, . . . , 6 and g n = 0 otherwise. Genes of group a are defined as showing no effect on the EGF treatment and therefore g n equals 0 by definition.
By definition, we obtain that g 1 = 0 for each of the three groupsz. From the first model assumption we obtain that g 1 , g 3 , and g 5 are equal to 0 for each of the three groupsz. From the second model assumption we obtain that (g 2 , g 4 , g 6 ) is equal to the corresponding column of Fig. 5 for each of the three groupsz. Figure 6 summarizes the values of the indicator variables g 1 , . . . , g 6 for each of the three groups b − d.
Third, we assume that the logarithmic expression levels x 1 , . . . , x 6 are statistically independent and normally distributed. By combining all three model assumptions, we obtained the likelihood p(x|a, θ a ) = denotes the density of the normal distribution, θz = (μz 0 , μz 1 , σz) denotes the parameter of modelz, and g n are the indicator variables from Fig. 6.

Posterior approximation by the Bayesian Information Criterion
Next, we seek the approximate posterior for each z ∈ {a, b, c, d} and each gene, where p(z) is the prior probability of group z.
For the four models of section "Probabilistic modeling of gene expression" the approximations of the marginal