Comment on "The Molecular Evolutionary Patterns of the Insulin/FOXO Signaling Pathway".

Letter to the Editor on Wang M, Wang Q, Wang Z, Zhang X, Pan Y. The molecular evolutionary patterns of the insulin/FOXO signaling pathway. Evol Bioinform. 2013;9:1-16. doi: 10.4137/EBO.S10539.

Wang et al have recently analyzed the pattern of selective pressures acting on genes along the insulin/FOXO pathway of vertebrates.1 In this study, which is similar to those that we and other authors performed in the Insulin/TOR (IT) signal transduction pathway,2–5 the authors questioned some of our previous findings. Here we revise previous analyses of the pattern of selective constraints along the IT and other signal transduction pathways and discuss a number of issues that might account for the apparently different results obtained for the IT pathway by the different groups.

Alvarez-Ponce et al2 analyzed the selective pressures acting on genes involved in the Drosophila IT pathway and found that levels of selective constraint exhibit a gradient distribution along the upstream/downstream pathway axis. In particular, levels of selective constraint (measured from the nonsynonymous to synonymous divergence ratio, ω = dN/dS) correlate with pathway position (computed as the number of steps required to transduce the signal from the insulin receptor to the rest of pathway components), being the downstream genes the most selectively constrained. We also determined that this trend is not the result of a higher incidence of positive selection in the upstream part of the pathway2,4 and that the correlation remains significant after factoring out the effect of gene expression level, codon bias, protein length, and number of protein-protein interactions. The analysis was based on 27 D. melanogaster genes (representative of 21 gene families), all of which act downstream of the insulin receptor (ie, neither insulin-like peptides, IGFBPs, nor the insulin receptor were included in this analysis). We identified the orthologs of each gene in another 11 Drosophila species and, when necessary, we corrected manually the gene annotations. Remarkably, the distribution of selective constraint levels observed along the IT pathway sharply contrasts with that generally observed in biosynthetic pathways, in which upstream genes tend to be the most selectively constrained.6–11

In order to evaluate whether the correlation observed in Drosophila is general or, on the contrary, specific to this genus, we performed a similar analysis using six vertebrate genomes.3 In this case, we analyzed a total of 72 genes (representative of 23 gene families) (Table 1). Our network-level analysis revealed a tendency consistent with that observed in Drosophila: levels of selective constraint exhibited a polarity along the pathway, with downstream genes being the most constrained. This trend, however, was less pronounced than that observed in Drosophila, which might be explained, at least partially, by the lower effective population sizes that vertebrates commonly have.12,13 In fact, overall ω values (estimated using the M0 model implemented in the PAML software,14 which assumes a unique ω value for all codons and branches in the phylogeny) do not significantly correlate with pathway position; however, when the analysis is conducted separately for each branch of the phylogeny of the studied vertebrates (using the ω values estimated under the free-ratio model, which assumes a separate ω ratio for each branch), the sign of the correlation was negative (ρ < 0) for most of the branches. This result indicates that, also in vertebrates, levels of selective constraint follow a polarity along the pathway, thus implying that the trend does not represent a feature specific to Drosophila. Jovelin and Phillips also studied the levels and patterns of selective constraint along the IT pathway, in this case in Caenorhabditis.5 In their similar network-level analysis based on 13 genes, they also found that ω values significantly correlated with pathway position, with downstream genes being the most selectively constrained. Nevertheless, they observed that the trend becomes nonsignificant after discounting the effect of gene expression levels.

Table 1

Genes analyzed by Alvarez-Ponce et al,2,3 Jovelin and Phillips,5 and Wang et al.1

Gene family	Pathway positiona	Alvarez-Ponce et al. 2009 (Drosophila)2	Alvarez-Ponce et al. 2011 (vertebrates)3	Jovelin and Phillips 2011 (Caenorhabditis)5	Wang et al. 2013 (vertebrates)1
IGFBP	−1	–	–	–	6
INS/IGF1	−1	–	–	–	2
INSR	0	–	3	1	1
IRS	1	1	11	1	3
P85	2	1	3	1	3
P110	3	1	4	1	2
MELT	4	1	1	–	–
PDK1	4	1	2	1	–c
AKT	5	1	3	1	3
SGK	5	–	–	1	–
PKC	5	–	9	–	–
TSC1	5	1	1	–	–
FOXO	6	1	4	1	3
GSK3	6	1	2	–	–
TSC2	6	1	1	–	–
EIF2Bɛ	7	1	1	–	–
GYS	7	1	2	–	–
MYC	7	1	4	–	–
RHEB	7	1	2	–	–
TOR	8	1	1	1	–
RPTOR	8	–	–	1	–
4EBP	9	1	3	1	–
S6K	9	1	2	–	–
EIF4E	10	7	5	–	–
RPS6	10	1	2	–	–
FOXL1	10	–	–	1	–
HIF	10	–	–	1	–
CYH	–b	1	4	–	–
PTEN	–b	1	2	–	–
Total		27	72	13	23

Notes: Number of genes per family analyzed.

Number of steps required to transduce the signal from the insulin or IGF1 receptors to each protein in the pathway. The same numbering as in ref. 2 is used. In the other works, a consistent numbering is used;

gene families not included in network-level analyses, as they do not directly participate in the main signaling cascade;

Wang et al1 included in their analyses genes encoding pyruvate dehydrogenase kinases 1–4 and not the 3-phosphoinositide dependent protein kinase-1.

Noticeably, a similar distribution of selective constraint levels has been observed in the yeast HOG signal transduction pathway, whose downstream genes are the most selectively constrained.15 Indeed, this appears to be the most common distribution of levels of selective constraint across signal transduction cascades: analysis of the global human signal transduction network has shown that, on average, genes acting at the downstream part of the network tend to be more selectively constrained than those occupying more highly hierarchical positions.16

In their recent study in vertebrates, Wang et al1 used a gene set that only partially overlaps with that used in our 2011 analysis. Their gene set encompasses a total of 27 genes (representative of 9 gene families) acting at the upper part of the IT pathway (the insulin/FOXO pathway), of which only 15 (representative of 5 gene families) were also used in our study (Table 1).3 Consistent with previous observations, 2–5 they also detected a negative significant correlation between ω estimates and pathway position, even when potentially confounding factors were controlled for. However, they noticed that statistical significance vanished when genes encoding insulin, IGF1, and IGFBP1–6 (genes not included in any of the previous analyses2–5) were removed from their analysis. From this observation, the authors conclude that the distribution of selective constraints that they observe across the pathway “is different from [our] previous similar reports.” They claim that their observation that the detected trend vanishes after removing the above-mentioned genes “demonstrates that the relaxed selective constraint at position 1 explains the purifying selection polarity along the IIF pathway,” adding that “the presented correlation [between pathway position and ω] comes from the relaxed selective constraint operated on genes of position 1, which give us an erroneous perception of truth.”1 Wang et al’s conclusions, however, are not supported by their observations, as none of the previous analyses included genes encoding insulin, IGF1, and IGFBP1–6 (genes at position 1 in their study) (Table 1).2–5 Therefore, the polarity that we had detected cannot be caused by the presence of these genes.

Furthermore, we have reanalyzed our published data sets to evaluate whether the correlation between ω and the IT pathway position previously detected vanishes when analyses are restricted to the gene families also analyzed by Wang et al (Table 1). Remarkably, when we reanalyzed our Drosophila dataset, we found a significant correlation that is even stronger than that detected in our previous analysis of the entire pathway2 (Spearman rank correlation coefficient, ρ = −0.898, P = 0.015; Kendall rank correlation coefficient, τ = −0.828, P = 0.023) (Fig. 1). The reanalysis of our vertebrate dataset shows that the correlation between ω and pathway position is negative (negative Spearman and Kendall rank correlation coefficients) in 8 out of the 9 branches of the phylogeny (Table 2). This number of branches with a negative correlation coefficient is higher than the 50% expected under a random distribution of levels of selective constraint along the pathway (binomial test, P = 0.039). Therefore, current results in Drosophila and vertebrates are consistent with our previously published results and indicate that our previous findings were not the result of using a particular set of genes, since the trend can also be found in the specific part of the pathway studied by Wang et al.1

Figure 1

Correlation between ω estimates and pathway position in the Drosophila insulin/TOR pathway. The analysis is restricted to the gene families used both by Alvarez-Ponce et al2 and by Wang et al.1 The ω estimates were obtained from Alvarez-Ponce et al.2 Positions 1–6 correspond, respectively, to genes chico, Pi3K21B, Pi3K92E, Pk61C, Akt1 and foxo.

Table 2

Correlation between ω estimates and pathway position in the vertebrate insulin/TOR pathway.

Phylogenetic branch	ρ
a (Homo sapiens)	0.267
b (Mus musculus)	−0.272
c (Bos taurus)	−0.195
d (Monodelphis domestica)	−0.208
e (Ornithorhynchus anatinus)	−0.011
f (Gallus gallus)	−0.362
g (internal branch)	−0.336
h (internal branch)	−0.348
i (internal branch)	−0.497

Notes: Analysis restricted to the gene families used both by Alvarez-Ponce et al3 and by Wang et al.1 Estimates of ω were obtained from Alvarez-Ponce et al (Dataset 1).3 Branches are named as in Alvarez-Ponce et al.3

Here we highlight the most important differences between the studies conducted by the different groups, which might account for the observed dissimilarities: quality of gene annotations, genes analyzed, and genes not involved in the main insulin signaling cascade.

Quality of Gene Annotations

In our previous analyses, we visually inspected all multiple sequence alignments to identify putative gene annotation errors (eg, unannotated exons or exons erroneously annotated in some species). Such annotation errors were, wherever possible, manually fixed using the information of the original genome sequences. After this manual improvement of gene annotations, alignment regions that remained problematic were removed from the analysis. Jovelin and Phillips5 also manually annotated exon/intron structures.5 Although Wang et al manually removed “unreliably aligned regions,”1 they do not mention having corrected annotation errors through manual annotation. Because genome annotations of non-model organisms are mostly the output of automatic computational pipelines, they generally contain a high proportion of errors (See, for example, Tu et al17 and Collins et al.).18 Manual correction of such errors is very important, since ω estimates are highly sensitive to errors in the gene structure annotation (See, for example, Schneider et al.).19

Genes Analyzed

Wang et al included only 27 genes in their analysis (representative of 9 gene families), all of which act at the upper part of the vertebrate IT pathway, whereas in our vertebrate analysis we used a total of 72 genes (belonging to 23 families), with only 15 (5 families) being shared between both datasets (Table 1).1,3 Their dataset includes genes encoding insulin, insulin-like growth factor 1 (IGF1), and IGF-binding proteins 1–6 (IGFBP1–6), which were not included in our previous analyses. Conversely, their dataset does not contain genes encoding proteins IGF1R, DOK, FRS, p110δ, PKC, Melted, GSK3, TSC1/2, Rheb, mTOR, S6K, RpS6, 4E-BP, eIF4E, Myc, eIF2B-ɛ, Glycogen Synthase, PTEN, and Cyh, which were included in our previous analysis (Table 1). The fact that they used a reduced dataset may have decreased the statistical power of their correlation analyses, particularly when genes encoding insulin, IGF1 and IGFBP1–6 (ie, 30% of their dataset) were removed from their analyses. Notice that trimming a small dataset often involves a dramatic reduction in the power of statistical tests, which can result in certain trends becoming statistically nonsignificant, despite their relevance.

Genes not Involved in the main Insulin Signaling Cascade

According to their report, Wang et al1 included in their analyses 4 genes encoding 3-phosphoinositide dependent protein kinases (“position 5” in their Fig. 1, downstream of the PI3K complex and upstream of AKT1-3 proteins). However, the accession numbers listed in their additional files 1 and 2 correspond to genes encoding pyruvate dehydrogenase kinases 1–4. Although the activity of the latter enzymes is regulated by insulin, they do not act at “position 5” (in their Fig. 1); indeed, they are not directly involved in the main insulin signal transduction cascade. While insulin activates the 3-phosphoinositide dependent protein kinase-1, it inhibits pyruvate dehydrogenase kinases 1–4.20–23 The fact that the same acronym (PDK) is often used for both enzymes (even though the gene symbols are PDPK1 for the 3-phosphoinositide dependent protein kinase-1 and PDK1–4 for pyruvate dehydrogenase kinases 1–4) likely led Wang et al to confound both gene families.

Conclusions

In conclusion, the issues raised above might explain the observation by Wang et al1 that the negative correlation between pathway position and ω vanishes after removing some genes from the analysis. These concerns, indeed, may have added some noise to their dataset, reducing the statistical power of their analyses. Therefore, there is no support for claiming that levels of selective constraint do not exhibit a polarity along the upstream/downstream axis of the insulin/TOR pathway.