s2fieldsofstudy: Biology Medicine

text:
Thioredoxin interacting protein (TXNIP) is a novel tumor suppressor in thyroid cancer

Background Thyroid cancer is the most common endocrine malignancy, and many patients with metastatic differentiated thyroid cancer (DTC), poorly differentiated thyroid cancer (PDTC), and anaplastic thyroid cancer (ATC) fail to respond to conventional therapies, resulting in morbidity and mortality. Additional therapeutic targets and treatment options are needed for these patients. We recently reported that peroxisome proliferator-activated receptor gamma (PPARγ) is highly expressed in ATC and confers an aggressive phenotype when overexpressed in DTC cells. Methods Microarray analysis was used to identify downstream targets of PPARγ in ATC cells. Western blot analysis and immunohistochemistry (IHC) were used to assess thioredoxin interacting protein (TXNIP) expression in thyroid cancer cell lines and primary tumor specimens. Retroviral transduction was used to generate ATC cell lines that overexpress TXNIP, and assays that assess glucose uptake, viable cell proliferation, and invasion were used to characterize the in vitro properties of these cells. An orthotopic thyroid cancer mouse model was used to assess the effect of TXNIP overexpression in ATC cell lines in vivo. Results Using microarray analysis, we show that TXNIP is highly upregulated when PPARγ is depleted from ATC cells. Using Western blot analysis and IHC, we show that DTC and ATC cells exhibit differential TXNIP expression patterns. DTC cell lines and patient tumors have high TXNIP expression in contrast to low or absent expression in ATC cell lines and tumors. Overexpression of TXNIP decreases the growth of HTh74 cells compared to vector controls and inhibits glucose uptake in the ATC cell lines HTh74 and T238. Importantly, TXNIP overexpression in T238 cells results in attenuated tumor growth and decreased metastasis in an orthotopic thyroid cancer mouse model. Conclusions Our findings indicate that TXNIP functions as a tumor suppressor in thyroid cells, and its downregulation is likely important in the transition from differentiated to advanced thyroid cancer. These studies underscore the potential of TXNIP as a novel therapeutic target and prognostic indicator in advanced thyroid cancer.

Background
Thyroid cancer is the most common endocrine malignancy, and it is estimated that nearly 63,000 new cases of thyroid cancer will be diagnosed in the United States in 2014 [1]. Although the majority of these patients have welldifferentiated thyroid cancer and respond favorably to conventional therapies (surgery with or without radioactive iodine I 131 therapy and suppression therapy with thyroid hormone), a significant minority of patients develops advanced disease that is resistant to standard treatments. For the minority of individuals who develop anaplastic thyroid cancer (ATC), an aggressive and undifferentiated form of thyroid cancer, the prognosis is very poor and median survival is 3-5 months [2,3]. The mechanisms underlying the development of poorly differentiated thyroid cancer (PDTC) and ATC are incompletely understood. Therapeutic options for these patients are limited, and prognosis remains dismal. This is an area for which further research and drug/therapy development is critically needed.
We recently reported that peroxisome proliferatoractivated receptor gamma (PPARγ) confers an aggressive phenotype in thyroid cancer cells [4]. Nuclear PPARγ expression was absent in differentiated thyroid cancer (DTC) cell lines and high in ATC cell lines. When PPARγ was overexpressed in the DTC cell line BCPAP, in vitro proliferation and invasive capacity were increased. Furthermore, when PPARγ was depleted from ATC cells, in vitro proliferation and invasive capacity were inhibited and tumor growth was inhibited in two in vivo murine cancer models (an orthotopic thyroid model and a flank xenograft model). These data challenge the widelyheld assumption that PPARγ is a tumor suppressor and suggest that PPARγ may mediate the aggressive phenotype that develops as part of the transition from DTC to PDTC and ATC.
In the studies presented here, we show that TXNIP, the gene encoding thioredoxin interacting protein (TXNIP), is a negatively-regulated downstream target of PPARγ. TXNIP is a negative regulator of cell growth and metabolism [5][6][7][8][9][10][11]. It modulates cellular redox status by binding to and inhibiting thioredoxin, a principal component of the cell's antioxidant system [12][13][14]. Furthermore, via its negative regulation of thioredoxin, TXNIP can inhibit invasion and metastasis and promote a pro-apoptotic cellular environment [15][16][17][18][19][20]. TXNIP has been shown to be a tumor suppressor in cancer [8,[21][22][23][24][25][26][27][28], but its role in thyroid cells or in thyroid cancer has not been investigated. In this current study, we show that TXNIP is highly expressed in DTC and low or undetectable in ATC and appears to be a novel tumor suppressor in thyroid cancer.

TXNIP is upregulated in PPARγ-depleted ATC cells
We previously reported that PPARγ confers an aggressive phenotype in thyroid cancer cells [4]. To investigate downstream mediators of these effects, a discovery-based microarray approach was used to compare mRNA expression levels in HTh74 ATC cells expressing a PPARγ-specific shRNA compared to a scrambled control. The most highly up-and down-regulated genes are shown in Table 1. PPARG was the most highly downregulated gene, which is consistent with the use of a PPARγ-specific shRNA. The most highly upregulated gene in the PPARγdepleted HTh74 cells was thyroid peroxidase, an enzyme critical to synthesis of thyroid hormone. Interestingly, TXNIP was the second most highly upregulated gene in the PPARγ-depleted ATC cells (>10-fold). TXNIP is a known tumor suppressor whose role in thyroid cancer has never been reported. Western blot analysis of whole cell extracts confirms that TXNIP is also upregulated at the protein level when PPARγ expression is knocked down in the HTh74 ATC cells ( Figure 1A). These data imply that TXNIP downregulation in HTh74 cells is a downstream consequence of high nuclear PPARγ expression, as TXNIP expression is increased when HTh74 cells are depleted of PPARγ. These data are consistent with a previously published report that the TXNIP promoter contains multiple PPARγ binding sites and that TXNIP expression is negatively regulated by binding of PPARγ to the TXNIP promoter and by treatment with PPARγ agonists [29].
TXNIP is expressed at high levels in DTC cells and tissues compared with ATC, and TXNIP expression correlates inversely with glucose uptake We next investigated whether there is differential TXNIP expression in DTC versus ATC cells. Western blot analysis revealed high TXNIP expression in the DTC cell lines TPC-1, MDA-T41, and BCPAP in contrast to low or undetectable expression in the ATC cell lines HTh74, C643, Ocut-2, and 8505C ( Figure 1B). This disparate TXNIP expression between DTC and ATC cell lines, however, was not absolute. The ATC cell lines T238 and TJH11T had levels of TXNIP expression comparable to some DTC cells, whereas the DTC cell line K1 had lower TXNIP expression in comparison to the other DTC cell lines. Similar trends were observed at the mRNA level (data not shown).
TXNIP is a known binding partner and inhibitor of the cellular antioxidant and tumor promoter thioredoxin (Trx) [12,14]. Thioredoxin may be localized to the cytoplasm or nucleus (Trx-1) where it can function as a cofactor, detoxify reactive oxygen species, or modulate the activity of transcription factors, or to mitochondria (Trx-2), and high Trx-1 expression is associated with more aggressive disease in other cancer types [30][31][32][33][34][35][36]. As TXNIP is a negative regulator of Trx and disparate TXNIP expression was observed between the panels of DTC and ATC cell lines, we sought to identify differences in Trx-1 expression patterns between the two groups. Immunoblot analysis revealed no clear differences in Trx-1 expression levels between DTC and ATC ( Figure 1B).
TXNIP is also a known inhibitor of glucose uptake [5,7,11]. We, therefore, investigated whether differential TXNIP expression correlates with differences in glucose uptake. Analysis of glucose uptake of two ATC and two DTC cell lines was performed, and the ATC cell line HTh74, which has the lowest TXNIP levels of the examined cell lines, had the highest level of glucose uptake ( Figure 1C). Interestingly, the ATC cell line T238 had glucose uptake at levels comparable to that observed with the two DTC cell lines TPC-1 and MDA-T41. This finding is consistent with a higher level of TXNIP expression in this cell line ( Figure 1B). We next investigated whether TXNIP was differentially expressed at the tumor tissue level. Immunohistochemistry to detect TXNIP protein expression was performed on paraffin-embedded tumor blocks of 13 well-differentiated primary papillary thyroid cancer (PTC) tumors and 8 ATC patient specimens. Representative images at high and low magnification are shown for PTC (Figure 2A-B) and ATC ( Figure 2C). TXNIP expression was undetectable in 63% of ATC specimens in contrast to 15% of PTC specimens, though there was one ATC specimen that exhibited high TXNIP staining. This differential TXNIP expression between PTC and ATC is similar to the pattern observed with the thyroid cancer cell lines ( Figure 1B). Therefore, loss of or absence of TXNIP expression appears to correlate with more aggressive thyroid cancer in cell lines and tumor tissue.

TXNIP overexpression in the ATC cells attenuates in vitro glucose uptake and growth
Based on its differential expression in DTC and ATC, we predicted that TXNIP acts as a tumor suppressor in thyroid cells and that its downregulation plays an important role in the development of an aggressive thyroid cancer phenotype. To further investigate this hypothesis, we reexpressed TXNIP in ATC cell lines. HTh74 and T238 ATC cell lines were transduced with a retroviral vector encoding human TXNIP. Western blot analysis verified increased TXNIP expression in the transduced HTh74 ( Figure 3A) and T238 ( Figure 3B) cell lines compared to control cells transduced with empty vector. As a functional read out of TXNIP expression, we assessed glucose uptake in the stable cell lines. TXNIP overexpression significantly inhibited glucose uptake in both the HTh74 and T238 cell lines ( Figure 3C and 3D, respectively), consistent with its known function of glucose uptake inhibition in other tissues.
To assess the effect of TXNIP overexpression on the growth of these ATC cell lines, viable cell proliferation assays were performed under standard growth conditions. TXNIP overexpression resulted in slowed growth of HTh74 cells by 37% ( Figure 4A). Interestingly, however, TXNIP overexpression in T238 cells had no effect on the in vitro growth rate ( Figure 4B). The reasons for the observed differences in TXNIP-mediated growth effects between these two cell lines are unclear, however, there are a few potential causes or factors that might contribute to this discrepancy. The baseline proliferation rate of parental T238 cells is much higher than parental HTh74 cells (data not shown), and this may obscure our ability to detect more subtle growth inhibitory effects with TXNIP overexpression in the T238 cell line. Alternatively, differential pathway activation might explain the in vitro growth differences between T238 and HTh74 cells and the ability of T238 cells to circumvent TXNIP-mediated growth deceleration. The T238 parental cell line has some basal TXNIP expression, and these cells have likely acquired the ability to resist some of the in vitro growth inhibitory effects of TXNIP via other mechanisms. In addition, the standard growth media used in propagation of these cell lines provides supplemental glutamine and glucose, and this enriched nutrient media may permit Stable cell lines were generated via transduction of the ATC cell line HTh74 with lentivirus encoding PPARγ-specific shRNA or scrambled control followed by selection in puromycin 0.5 μg/mL. Prepared RNA was then subjected to microarray analysis. Genes whose mRNA levels varied by greater than 5-fold are shown in the table along with fold change.
T238 cells to overcome the metabolic inhibitory effects of TXNIP overexpression in vitro. These potential limitations of in vitro assay systems, which do not adequately recapitulate in vivo tumor conditions, underscore the importance of further evaluation in an in vivo model system. In addition to proliferation assays, we performed invasion assays using an in vitro Matrigel invasion model [37], and TXNIP overexpression resulted in a trend towards decreased invasion in both cell lines but this did not reach statistical significance ( Figure 4C-D).

TXNIP expression in ATC cells results in attenuated tumor growth and metastasis in an in vivo orthotopic thyroid cancer mouse model
Finally, we examined the effect of TXNIP overexpression in an in vivo orthotopic tumor model. The orthotopic thyroid cancer mouse model is a well-established model that closely mimics the features of human thyroid cancer with regard to growth and metastases than does the more commonlyused subcutaneous flank xenograft model [4,[38][39][40][41]. Tumor cells expressing luciferase-IRES-GFP were injected into the right thyroid lobe and monitored weekly by IVIS imaging for tumor establishment and growth. TXNIP overexpression in the T238 cell line resulted in attenuated bioluminescence compared to vector control ( Figure 5A-B), and resultant tumor volumes were significantly smaller ( Figure 5C). Though there was a significant attenuation in bioluminescence as well as a trend toward smaller tumors with the TXNIP-expressing HTh74 cells compared to vector controls, the final tumor volumes were not significantly different (data not shown). In our experience, T238 cells in this orthotopic model frequently result in lung metastases (unpublished data). To determine if TXNIP overexpression affected distant metastatic spread in this orthotopic model, lungs were collected at necropsy and snap-frozen, homogenized, and RNA was subsequently isolated and subjected to qRT-PCR to detect GFP expression, a marker of lung metastases in this model. Metastatic tumor burden was significantly reduced in the mice injected with TXNIP- overexpressing T238 cells compared to vector controls ( Figure 5D). Control mice with no tumor cells injected had no detectable GFP mRNA in their lung tissues.

Discussion
In this report, we identified TXNIP as a novel tumor suppressor in thyroid cancer. TXNIP is highly upregulated when PPARγ, which we have previously shown to be a tumor promoter [4], is depleted from ATC cells. Furthermore, DTC cell lines and primary PTC tumors have high endogenous TXNIP levels, whereas TXNIP expression is low or absent in ATC cell lines and primary tumor specimens. This TXNIP expression pattern is opposite to what we previously reported with PPARγ, which is highly expressed in ATC cell lines and absent in DTC cell lines and whose forced expression confers a more aggressive phenotype in thyroid cancer cells in vitro and in vivo [4]. These data are consistent with a previously published report that TXNIP is a negatively-regulated target of PPARγ [29]. Therefore, TXNIP expression appears to be downregulated or lost in the progression from well-differentiated thyroid tumors to more aggressive, undifferentiated tumors. TXNIP has been shown to be a tumor suppressor in other cell types, and our results here show for the first time that it serves a similar function in thyroid cells.
The novel finding that TXNIP expression is lost in the progression from well-differentiated PTC to undifferentiated, aggressive ATC is consistent with our hypothesis that TXNIP is a tumor suppressor in thyroid cancer. Patients with well-differentiated PTC respond well to conventional therapy and have an excellent overall survival rate. Poorly-differentiated PTC tumors that have lost the ability to concentrate iodine often fail to respond to conventional therapy and result in poorer outcomes, and undifferentiated ATC portends an extremely-poor prognosis and is generally fatal within 3-5 months [2,3]. The apparent loss of TXNIP expression during the progression from well-differentiated to poorly-differentiated and undifferentiated thyroid cancer is consistent with its role as a tumor suppressor in thyroid cells. Loss of TXNIP expression has been reported to correlate with more aggressive disease, advanced stage, and poorer prognosis in breast, gastric, colorectal, and bladder cancers, as well as diffuse large B-cell lymphoma [23,24,[42][43][44][45][46][47], and TXNIP mRNA expression has been shown to be inversely proportional to melanoma progression [27]. TXNIP expression in cancer may be downregulated through epigenetic, transcriptional, post-transcriptional, or translational mechanisms (reviewed by Zhou et al. [21]). Though the tumor promoter and cellular antioxidant Trx-1, which is inhibited by TXNIP, has been shown to confer more aggressive disease in other cancers [30][31][32][33][34], our data failed to show differences in Trx-1 expression levels between DTC and ATC cells.
Overexpression of TXNIP in the HTh74 ATC cell line resulted in slowed in vitro growth. This negative growth regulatory effect has been seen in other systems as well. TXNIP overexpression in the human gastric carcinoma cell lines AGS, SNU-16, and SNU-620, the promyelocytic leukemia cell line HL-60, and HTLV-I-positive T cells led to growth reduction in vitro [8,22,48]. Lung fibroblasts from TXNIP knockout mice proliferate at a faster rate than wild-type, implying that loss of TXNIP promotes or allows for enhanced proliferation [9]. Interestingly, we did not observe an in vitro growth inhibitory effect in the T238 cell line, which has some basal endogenous TXNIP expression. It is likely that the effects of TXNIP on cell growth and proliferation are cell-context dependent and might be circumvented through activation of alternative mitogenic pathways. Furthermore, in vitro cell culture conditions do not adequately recapitulate the tumor microenvironment and contributions of paracrinemediated signaling, underscoring the importance of in vivo studies. In accordance with this potential limitation, Goldberg and colleagues failed to see slowed in vitro growth of melanoma cells transfected with TXNIP though when injected in an orthotopic flank model in nude mice, slowed tumor growth/development was observed [27]. Although we observed a trend toward decreased invasion by TXNIP overexpression in two ATC cell lines, this effect did not reach statistical significance in our in vitro model. Figure 2 TXNIP is expressed in primary DTC tumors but is low or undetectable in ATC tumors. Immunohistochemistry was performed using a TXNIP-specific primary antibody or normal mouse IgG control and a horseradish peroxidase tagged secondary antibody on a panel of 13 PTC and 8 ATC patient specimens. Two representative PTC (A-B) and one ATC (C) specimen stains are shown at 20X and 60X power. Positive staining is indicated by brown coloring, though staining of colloid is nonspecific. Staining was quantitated by a pathologist (S. B. S.), and intensity of staining was scored from absent (0) to high to (3+), followed by percent of tumor positivity.
Importantly, in a well-established orthotopic murine thyroid cancer model that mimics human thyroid cancer with regard to growth and metastasis, we show that TXNIP overexpression in the ATC T238 cell line resulted in significant attenuation of both tumor growth and pulmonary metastatic burden. These data support our hypothesis that TXNIP is a tumor suppressor in thyroid cells. TXNIP has been shown to be a tumor suppressor in other animal models of cancer as well. A mouse strain with a spontaneous nonsense mutation in TXNIP has dramatically increased incidence of spontaneous hepatocellular carcinomas (HCC) [25], and TXNIP-knockout mice develop increased number and size of HCC in a diethylnitrosamine (DEN)-induced murine model of HCC [28]. In a murine gastric carcinoma model in which tumors are induced via infection with Helicobacter pylori and cotreatment with N-methyl-N-nitrosourea, concomitant knock out of TXNIP resulted in increased numbers of tumors, heightened preneoplastic changes, increased percentage of malignant tumors, and elevated inflammatory marker expression compared to control mice with wild-type TXNIP expression [22]. In a murine model of bladder cancer in which tumors are induced by treatment with N-butyl-N-(4-hydroxybutyl) nitrosamine (BBN), genetic deletion of TXNIP results in accelerated development of high grade and invasive tumors by~4 weeks compared to controls with wild-type expression, however, controls eventually succumb to tumor development and TXNIP expression in these tumors has been downregulated by other mechanisms [24].
In the orthotopic ATC model, TXNIP overexpression also led to a significant reduction in pulmonary metastatic burden. Inhibition of metastasis conferred by TXNIP overexpression has been shown in other systems as well. B16F10 melanoma cells transfected with TXNIP then injected via tail vein into C57BL/6 mice resulted in decreased lung metastases [23]. TXNIP-transfected melanoma cells resulted in fewer metastases in both a nude mouse flank tumor model and IV tail injection metastasis model relative to vector controls [27]. In human breast cancer, high TXNIP levels are associated with longer metastasis-free intervals and better prognosis than those with low TXNIP expression [43,46]. These data implicate TXNIP as a tumor suppressor in a variety of cancers and, for the first time, is now shown to be a tumor suppressor in thyroid cells.
Curiously, TXNIP overexpression in the ATC cell line HTh74 resulted in reduced in vitro growth but no significant difference on in vivo growth in the orthotopic thyroid cancer model. Although bioluminescence signals were attenuated in the TXNIP-expressing HTh74 cells versus controls, final tumor volumes were not significantly Glucose uptake assays were performed as described in Figure 1 using the HTh74 stable cell lines (C) and T238 stable cell lines (D). Data from all experiments were combined, and glucose uptake from each cell line was normalized to vector control levels (average set at 1). **p = 0.001, ***p <0.0001. different, though a trend toward smaller tumors with injection of the TXNIP-expressing HTh74 cells compared to vector controls was observed. TXNIP-overexpressing and vector control tumors did not look different histologically. In our prior studies using HTh74 cells in the orthotopic murine thyroid cancer model system, we have observed that the in vivo growth rates are slower compared to other thyroid cancer cell lines (84 days to achieve 100 mm 3 tumors compared with 28-35 days in other ATC cell lines). It is possible that if our study had been temporally extended, the trend in tumor volume attenuation in the TXNIP-overexpressing group might have reached statistical significance. It is also possible that TXNIP does not play a significant in vivo role on malignant behavior in the HTh74 cells, but our in vitro data would suggest otherwise.
In keeping with the known function of TXNIP as a glucose uptake inhibitor, we showed that the degree of glucose uptake was inversely correlated with TXNIP levels in the examined thyroid cancer cell lines. An interesting aspect of thyroid cancer biology relates to its properties on 2-deoxy-2-fluoro-D-glucose positron emission topography computed topography (FDG PET/CT) imaging. FDG PET imaging can be negative in many patients with DTC and distant metastases, and this correlates with a relatively good prognosis in these patients [49,50]. Less differentiated PTC and ATC tumors are more likely to be PET positive, PET positive lesions are more likely to be resistant to conventional radioactive I 131 treatment, and increased intensity of FDG uptake is associated with a poorer prognosis and increased mortality [49,50]. The mechanism underlying this differential glucose uptake between well-differentiated PTC and ATC is not well understood. The novel finding that TXNIP expression is low in ATC is consistent with the observed FDG uptake on PET/CT in patients with ATC, supporting a critical role for TXNIP as a metabolic regulator in thyroid cancer progression.
In addition to inducing a metabolic shift important to tumor biology, downregulation of TXNIP has other important effects in cancer cells that contribute to tumor promotion and/or progression. TXNIP can reduce tumor invasion and angiogenesis through inhibition of thioredoxin and can directly impact cell survival by promoting a pro-apoptotic environment [13,[15][16][17][18][19][20]. Independent of its interaction with thioredoxin, TXNIP also has the ability to inhibit cell cycle progression by indirectly stabilizing the cell cycle inhibitor p27 Kip1 [9]. In addition, TXNIP indirectly inhibits mTOR, a regulator of cell growth and metabolism [6]. Therefore, downregulation of TXNIP in a tumor has the potential to promote cell survival, growth, invasion, and metastasis. The exact mechanisms by which TXNIP exerts its tumor suppressive functions in thyroid cancer cells are not yet clear. Future studies of the mechanisms by which TXNIP is expressed and functions in thyroid cancer will improve our understanding of the progression to advanced thyroid cancer and help to develop more effective targeted therapies.

Conclusions
In conclusion, we report that TXNIP is a novel tumor suppressor in thyroid cancer. TXNIP is downregulated during the progression from well-differentiated thyroid cancers to poorly differentiated and anaplastic thyroid cancers. Overexpression of TXNIP in ATC cell lines resulted in slowed in vitro growth and glucose uptake inhibition. Importantly, in an in vivo orthotopic murine thyroid cancer model, TXNIP overexpression attenuated tumor growth and drastically diminished pulmonary metastatic tumor burden. These data highlight the importance of TXNIP as a potential therapeutic target and prognostic marker in advanced thyroid cancer.

Cell lines and maintenance
HTh74 and C643 cells were obtained from Dr. K. Ain (University of Kentucky, Lexington, KY) with permission from N. E. Heldin (University Hospital, Uppsala, Sweden). To assess for lung metastasis, RNA isolated from whole lungs was subjected to qRT-PCR analysis to assess GFP expression. Data was normalized to 18 s RNA levels, and averages plus SD were plotted (D). RNA isolated from lungs of mice receiving no tumor cell injections was included as an additional negative control for GFP expression and is indicated with closed circles (•). *p <0.05, **p =0.0021, ***p <0.001, ****p <0.0001. passaged at 37°C in 5% CO 2 . Cell lines were authenticated by short tandem repeat (STR) profiling as previously described [51].

PPARγ knockdown and microarray analysis
PPARγ-depleted HTh74 cells and scrambled control cells were generated as previously described using lentivirus expressing PPAR-specific shRNA or scrambled control [4]. Total RNA from PPARγ-depleted and scrambled control cells was isolated using an RNeasy Mini Kit (Qiagen) according to the manufacturer's instructions. Integrity of the RNA preparation was verified on an Agilent Bioanalyzer 2100. Total RNA (5 μg) of each cell line was used for microarray analysis using the Human Genome U133 Plus 2.0 Array (Affymetrix), performed by the Gene Expression Core of the University of Colorado Denver, Anschutz Medical Campus (Aurora, CO). Gene expression profiles were normalized by robust multichip analysis (RMA), differentially expressed genes were analyzed by fold-change, using a cut-off of 2-fold, 122 and 198 genes were found to be up and down-regulated in the knockdown line. Enrichment analysis of the gene list was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID) analysis software.

Immunohistochemistry
We retrospectively selected formalin-fixed, paraffinembedded blocks of primary 13 PTC and 8 ATC specimens from the University of Colorado Hospital pathology archives for analysis of TXNIP protein expression by immunohistochemistry. Institutional review board approval was obtained. Sections were deparaffinized in Histoclear, rehydrated, and antigen retrieval in 10 mM sodium citrate buffer with 0.05% Tween 20, pH 6.0, was performed in a Biocare Medical decloaking chamber at 120°C for 5 minutes. Endogenous peroxidase activity was quenched by incubation in 3% hydrogen peroxide for 30 minutes at room temperature. Tissues were blocked with 5% goat serum in phosphate-buffered saline (PBS) with 1% bovine serum albumin for 1 hour at room temperature. Slides were incubated overnight at 4°C in primary anti-TXNIP antibody at a concentration of 1:400 (mouse monoclonal antibody IgG1, clone JY2, MBL International, catalog # K0205) diluted in antibody dilution buffer (0.01 M PBS, pH 7.2 with 0.05% sodium azide) or normal mouse IgG at an equivalent concentration as negative control (Santa Cruz Biotechnology, catalog # sc-2025). Slides were incubated in secondary goat anti-mouse antibody conjugated to horseradish peroxidase diluted in PBS at a concentration of 1:400 (Dako, catalog # P0447) for 1 hour at room temperature. For visualization, slides were incubated at room temperature for 2-4 minutes in ImmPACT 3, 3′-diaminobenzidine (DAB) peroxidase substrate (Vector Laboratories). Sections were counterstained in Mayer's hematoxylin solution, dehydrated, dried, and mounted with Cytoseal (Thermo Scientific). Stained specimens were read and score by a pathologist (S. B. Sams) based on percent of specimen that stained positively and degree of intensity (0 to 3+).

TXNIP overexpression
A plasmid encoding human TXNIP (pcDNA3.1-hTXNIP) was a kind gift from P. Patwari (Brigham and Women's Hospital, Boston, MA). The TXNIP coding sequence was amplified via polymerase chain reaction (PCR), using primers TXNIP 1 F [5′ TAG CGG CCG CAT GGT GAT GTT CAA GAA GAT CAA GT 3′] and hTXNIP EcoRI rev [5′ GCG AAT TCT CAC TGC ACA TTG TTG TTG AGG A 3′], which added a NotI restriction site to the 5′ end and an EcoRI site on the 3′ end of the coding sequence, respectively. PCR product was ligated into the pCR2.1 shuttling vector (Invitrogen), excised via NotI and EcoRI digestion, gel purified, and directionally inserted into the retroviral vector pQCXIP (gift of S. Nordeen, University of Colorado Denver, Aurora, CO) pre-digested with NotI and EcoRI. Resultant pQCXIP-hTXNIP clones were confirmed with sequencing. Retrovirus was generated in the BOSC cell line [52] after transfection of pQCXIP or pQCXIP-hTXNIP with pCL-Ampho packaging vector (gift of H. Ford, University of Colorado Denver, Aurora, CO) using FuGENE 6 transfection reagent (Roche). Supernatant with virus was collected at 48 and 72 hours after transfection, centrifuged at low speed, filtered through 0.45 μm syringe filter (Fisher Scientific), and stored at −80°C. Anaplastic HTh74 and T238 cells were transduced with virus-containing supernatant mixed 1:1 with growth media and supplemented with 8 μg/mL polybrene (Sigma-Aldrich), as previously described [4]. In the vector pQCXIP, the coding sequence for the insert (TXNIP) is cotranscribed with a puromycin resistance gene as a bicistronic message via an internal ribosome entry site. Forty-eight hours after transduction, the cells were placed under selection in puromycin (Sigma-Aldrich) at a concentration of 0.5 μg/mL for HTh74 cells and 2.5 μg/mL for T238 cells as previously determined by kill curves.

Glucose uptake assays
Cells were grown in 12-well plates with each condition plated in triplicate. Prior to glucose uptake determination, cells were rinsed in PBS, then incubated in low glucose DMEM without serum for 4 hours at 37°C. Cells were then incubated in Krebs buffer (140 mM NaCl, 5 mM potassium chloride, 2.5 mM magnesium sulfate, 1 mM calcium chloride, 20 mM HEPES, pH 7.4) supplemented with dimethyl sulfoxide (DMSO) or 20 μM cytochalasin B (Sigma-Adrich), an actin polymerization inhibitor that blocks nonspecific glucose uptake, for 1 hour at 37°C. Next, the cells were incubated in 0.01 mM 2-deoxy-D-glucose [Sigma-Aldrich], 0.665 nCi/mL [1,2-3 H] 2-deoxy-D-glucose [PerkinElmer], and either DMSO or cytochalasin B (20 μM) in Krebs buffer for an additional 20 minutes at 37°C. After this time period, cells were immediately rinsed 3 times with ice-cold PBS, then lysed in 0.4 N sodium hydroxide. Base was subsequently neutralized with 0.4 N hydrochloric acid. Uptake of [ 3 H]2-deoxy-D-glucose was determined by scintillation counting (Beckman Coulter). Nonspecific glucose uptake as determined by the cytochalasin B group was subtracted, and glucose uptake in pmol was normalized to protein content as determined by the BioRad DC protein assay system. Experiments were performed at least 3 times with each cell line and condition in triplicate, and data were graphed and analyzed by t-test using GraphPad Prism software.

Viable cell proliferation assays
HTh74 and T238 cells stably expressing pQCXIP vector with or without TXNIP were plated in duplicate in 6 cm plates at 50,000 cells/plate in RPMI 1640 supplemented with 5% FBS, without antibiotics. At days 3, 5, and 7, cells were rinsed in PBS, incubated in 0.25% trypsin-EDTA, collected, and resuspended in RPMI with 5% FBS. Cells were counted via the ViCell automated cell counting system. On day 7, collected cells were subsequently lysed in EB and subjected to Western blot analysis to determine TXNIP protein expression. Experiments were performed at least 3 times, and data were combined, graphed, and analyzed by 2 way ANOVA using GraphPad Prism software.

Invasion assays
Invasion assays with 2×10 5 HTh74 and 1×10 5 T238 cells stably expressing QCXIP vector with or without TXNIP were performed as previously described using BD Biocoat Matrigel invasion chambers (8 μM pore size, 24-well; BD Biosciences) [4]. Five fields per well were counted using Metamorph software (Molecular Devices), and each condition was performed in triplicate. Data from three independent experiments were combined, and data averages were normalized to the vector control mean. Statistical analysis was performed via application of the two-tailed t-test using GraphPad Prism software.

Orthotopic tumor mouse model
The right thyroid lobes of athymic nude mice were injected with 500,000 T238 QCXIP and T238 TXNIP cells stably expressing a luciferase-IRES-GFP plasmid (pEGFP-Luc-N1, a kind gift from C. Li, University of Colorado Denver, Aurora, CO) in 5 μL PBS as previously described [4,[38][39][40][41]. Weekly bioluminescence imaging using Xenogen IVIS200 (Caliper Life Sciences) in the presence of injected luciferin substrate (Caliper Life Sciences) was performed to monitor tumor establishment and growth, and bioluminescence activity was analyzed using Living Image software (Xenogen Corporation). Bioluminescence curves were analyzed by 2-way ANOVA with Bonferroni posttests using GraphPad Prism software. There were 10-11 mice per group for each experiment, and the described experiment was performed two times. In toto, there were 21 mice in each experimental arm when data from the two independent studies were pooled. Animals were sacrificed at 26-28 days or sooner if ill or moribund, and final tumor dimensions were measured with calipers. Final tumor volumes were calculated using the formula (length × width × height)/0.5236 and compared with t-test using GraphPad Prism software. All procedures were conducted in accordance with a protocol approved by the Institutional Animal Care and Use Committee of the University of Colorado Denver.

Isolation of RNA from lungs and quantitative reverse transcription polymerase chain reaction (qRT-PCR) for eGFP expression
At time of sacrifice for the second mouse orthotopic injection experiment, lungs were collected, snap-frozen in liquid nitrogen, and stored at −80°C. Lungs from uninjected mice served as negative controls. To harvest RNA, lung tissues were diced in a petri dish on ice and homogenized in TRI Reagent (Sigma Aldrich) using sterile stainless steel beads and a Qiagen TissueLyser. Homogenized tissue in TRI Reagent (1 mL) was mixed with 200 μL chloroform, centrifuged, and aqueous phase (which contains RNA) was removed. RNA purification with Qiagen RNeasy kit was then performed per the manufacturer's instructions with the exception of an added step of column incubation with RNasefree DNase I stock solution (Qiagen) in between wash steps to remove any residual DNA. RNA was ultimately eluted with RNase free water, and RNA concentration was quantitated using a Synergy H1 microplate reader (BioTek).
GFP mRNA levels were measured by real-time qRT-PCR using an ABI Prism 7900 sequence detector (Applied Biosystems/Life Technologies). Primers and probe for GFP were designed with the assistance of the Prism 7900 sequence detection software (Primer Express, PE Applied Biosystems). The TaqMan probe was purchased from Life Technologies 5′ labeled with 6-carboxyfluo rescein (FAM) and 3′-labeled with 6-caboxy-tetramethyl rhodamine (TAMRA). The forward and reverse primer sequences were GFP-F 5′-CACATGGTCCTGCTGGA GTTC-3′ and GFP-R 5′-TTGTACAGCTCGTCCAT GCC-3′ and the TaqMan fluorogenic probe sequence was 6FAM-CCGCCGCCGGGATCACTCT-TAMRA. Amplification reactions were performed in MicroAmp optical plates (Applied Biosystems/Life Technologies) in a 20 μl mix containing 1X TaqMan Buffer A (500 mM KCl, 100 mM Tris-HCl, 0.1 M EDTA, 600 nM passive reference dye ROX, pH 8.3 at room temperature), 300 μM each of dATP, dGTP, dCTP and 600 μM dUTP, 5.5 mM magnesium chloride, 900 nM forward primer, 900 nM reverse primer, 200 nM probe, 1.25 U AmpliTaq Gold DNA Polymerase and the template cDNA. Thermal cycling conditions were as follows: 2 minutes at 50°C followed by activation of TaqGold at 95°C for 10 minutes. Subsequently, 40 cycles of amplification were performed at 95°C for 15 seconds and 60°C for 1 minute. Quantities of GFP in test samples were normalized to 18 s r-RNA (PE Applied Biosystems), and the Mann Whitney test was applied to assess for statistical significance using GraphPad Prism software.


extfieldsofstudy: Biology Medicine

text:
FISHER'S GEOMETRIC MODEL WITH A MOVING OPTIMUM

Fisher's geometric model has been widely used to study the effects of pleiotropy and organismic complexity on phenotypic adaptation. Here, we study a version of Fisher's model in which a population adapts to a gradually moving optimum. Key parameters are the rate of environmental change, the dimensionality of phenotype space, and the patterns of mutational and selectional correlations. We focus on the distribution of adaptive substitutions, that is, the multivariate distribution of the phenotypic effects of fixed beneficial mutations. Our main results are based on an “adaptive-walk approximation,” which is checked against individual-based simulations. We find that (1) the distribution of adaptive substitutions is strongly affected by the ecological dynamics and largely depends on a single composite parameter γ, which scales the rate of environmental change by the “adaptive potential” of the population; (2) the distribution of adaptive substitution reflects the shape of the fitness landscape if the environment changes slowly, whereas it mirrors the distribution of new mutations if the environment changes fast; (3) in contrast to classical models of adaptation assuming a constant optimum, with a moving optimum, more complex organisms evolve via larger adaptive steps.

Natural populations are constantly faced with environmental changes that force them to either adapt or go extinct. In Arabidopsis thaliana, Hancock et al. (2011) recently identified candidate single nucleotide polymorphisms scattered over the entire genome that affect flowering time and vernalization and are strongly correlated with climate variables. Likewise, annual cycles of reproduction of various plants and animals have been adjusted to the peak availability of food as a response to changing environments (Gienapp et al. 2013). Conversely, migratory bird species that fail to respond phenologically decline in population size (Møller et al. 2008). The brood parasitic common cuckoo (Cuculus canorus) population, for example, declined in size by 6% since 1980, as they failed to synchronize their reproductive and migratory cycles with those of their particular host species, to which they are adapted to in terms of egg size, coloration, and spottiness (Antonov et al. 2010;Møller et al. 2011).
In recent years, numerous theoretical studies of the population genetics of adaptation have attempted to provide a formal framework for the observed empirical phenomena (for a review see Orr 2005b). Central to these studies is the description of the fundamental event during adaptation, that is, the substitution of a resident allele (i.e., gene variant) by a beneficial mutation. The statistical description of this process has been at the heart of evolutionary biology (Charlesworth 1996), and is key to addressing seemingly simple questions, such as: From the set of mutations that emerge in a population, which are the ones that will get fixed and what is their effect on phenotype or fitness? Will adaptation proceed by many steps of small effect or just by a few adaptive substitutions of large effect? Do simple organisms evolve faster than complex ones?
One of the most influential models of adaptive phenotypic evolution is Fisher's geometric model (FGM) (Fisher 1930). In this model, a phenotype is treated as a point in a multidimensional trait space, and mutations are random vectors in this space, which are beneficial if they bring the mutant phenotype closer to a nearby local optimum. Thus, FGM implicitly assumes "universal pleiotropy" (each mutation affects every trait) and, therefore, equates pleiotropy with "organismic complexity." Despite its simplicity and the lack of a clear genetic context (Chevin et al. 2010), FGM, more than 80 years after its proposal, has yielded several robust predictions supported by growing empirical evidence: First, the distribution of fitness effects of new mutations is well approximated by a (displaced) negative gamma distribution (Martin and Lenormand 2006a; for empirical support see Hietpas et al. 2013). Second, the distribution of adaptive substitutions is approximately exponential, meaning that most fixed mutations are of small and only a few are of large effect (Orr 1998; for empirical support see Rockman 2012, but see Bell 2009). Finally, fixed mutational effects become on average smaller as organismic complexity increases (Orr 2000; for empirical support see Cooper et al. 2007)-a phenomenon that has been termed "the cost of complexity" (Orr 2000;Welch and Waxman 2003;Wagner and Zhang 2011).
The classical version of FGM, however, only addresses the situation in which a population is confronted with constant stabilizing selection after a sudden change in the environment (e.g., Orr 2002;Martin and Lenormand 2006a). In nature, in contrast, environmental change may as often be gradual (Hairston et al. 2005;Thompson 2005;Parmesan 2006;Perron et al. 2008). Collins (2011b) recently emphasized that "using [models of] instantaneous environmental change to understand adaptive evolutionary responses to gradual change will not only underestimate the amount of adaptation, but also predict the wrong genotypic and phenotypic changes." Indeed, the necessity to include gradual environmental change into studies of adaptive evolution has long been recognized in quantitative genetics (e.g., Maynard Smith 1976). A number of studies have focused on the so-called moving-optimum model, in which the optimal values of a quantitative trait change over time (Lynch and Lande 1993;Bürger and Lynch 1995;Waxman and Peck 1999;Bürger and Gimelfarb 2002;Nunney 2003;Collins et al. 2007;Gordo and Campos 2013); extensions include multivariate phenotypes and the effects of pleiotropic constraints (Jones et al. 2004;Gomulkiewicz and Houle 2009;Jones et al. 2012;Chevin 2013;Lourenço et al. 2013). The focus of these studies was, however, on the rate of adaptation (Lynch and Lande 1993;Bürger and Lynch 1995;Gomulkiewicz and Holt 1995;Nunney 2003;Hansen and Houle 2008;Chevin 2013;Kopp and Matuszewski 2014) and the evolution and maintenance of genetic variation (Bürger 1999;Waxman and Peck 1999;Bürger and Gimelfarb 2002;Jones et al. 2004Jones et al. , 2012Gomulkiewicz and Houle 2009). In contrast, characteristics of individual substitutions have been addressed only recently (Collins et al. 2007; Hermisson 2007, 2009a,b). In particular, Hermisson (2007, 2009a) employed the movingoptimum model of a single quantitative trait to study the fixation time of single mutations and the order in which mutations of dif-ferent phenotypic effect sizes become fixed. Their latest study (Kopp and Hermisson 2009b) addresses the distribution of adaptive substitutions during long-term adaptation. Specifically, they showed that this distribution is almost entirely determined by a scaled rate of environmental change γ, which combines ecological and genetic factors (see below), and is unimodal (with an intermediate mode) rather than exponential. That is, most substitutions have an intermediate phenotypic effect, while small-and large-effect substitutions are rare.
An obvious next question is how these results are affected if phenotypic adaptation to gradual change is constrained by pleiotropic correlations among the traits under selection (as frequently observed in nature; Svensson et al. 2001;Guerreiro et al. 2012;Roff and Fairbairn 2012). This is the aim of the present article. This way, we integrate two modeling traditions, which have had little overlap so far: on the one hand, the multivariate moving-optimum model as used by Jones et al. (2004Jones et al. ( , 2012, and on the other hand, Fisher's classical geometric model for the study of adaptive effect sizes (Fisher 1930;Orr 1998Orr , 2000. We study how the expected distribution of adaptive steps is influenced by the rate of environmental change, the number of traits under selection (i.e., "organismic complexity"), and by selectional and mutational correlations (i.e., the shapes of the fitness landscape and the multivariate distribution of new mutations). Our analysis shows that the genetic basis of adaptation crucially depends on the tempo and mode of environmental change.

Phenotype, environmental change, and selection
We consider the evolution of n phenotypic traits z = (z 1 , . . . , z n ) , each of which is under Gaussian stabilizing selection with regard to a time-dependent optimum z opt (t): where denotes transposition and (and thus also −1 ) is an n × n positive definite and symmetric matrix. Throughout this article, we choose the linearly moving optimum, where v = (v 1 , . . . , v n ) is the vector of environmental change. In the following, we will interchangeably refer to n as the "degree of pleiotropy" or the "degree of complexity." The matrix describes the shape of the fitness landscape (including a contribution of environmental noise to the phenotype z, which otherwise is not modeled explicitly; see Bürger 2000). We will say that selection is isotropic if is proportional to an identity matrix, = σ 2 I (σ 2 > 0); and selection is correlated if has nonzero off-diagonal entries. As a measure for the average width of the fitness landscape, we definē which is the geometric mean of the eigenvalues of (if the fitness landscape is represented by an ellipse, as in Figure S1_1 below, the axes of the ellipse have length proportional to the square root of the eigenvalues). Note that overall selection is strong ifσ 2 is small.

Genotypes and mutation
In accordance with Fisher's original model, we make the assumption of "universal pleiotropy," that is, each mutation affects every trait. We denote by α the vector of the phenotypic effects of a mutation, and we assume that its distribution p(α) (which we will refer to as the distribution of new mutations) is multivariate normal with mean 0 and covariance matrix M (thus, we assume a continuum-of-alleles model), that is, Like , M has dimensions n × n and must be symmetric and positive definite. The diagonal elements of M are the variances of the mutational effects for individual traits, whereas off-diagonal elements are the mutational covariances. We will say that mutation is isotropic if M is proportional to an identity matrix, M = m 2 I (m 2 > 0), and mutation is correlated if M has nonzero off-diagonal entries. A measure for the average variance of mutational effects (in an arbitrary direction) is given bȳ When comparing different degrees of pleiotropy/complexity, we typically assume that the distribution of mutational effects on a given trait is independent of the total number of traits n (so-called Euclidean superposition model ;Turelli 1985;Wagner 1988;Wagner and Zhang 2011). For example, with isotropic mutation (see above), adding more traits does not change the parameterm 2 . As a consequence, the average total effect of a mutation increases with n.
In Supporting Information 2, we introduce a transformation that shows that the general model outlined above can always be reduced to a model with isotropic selection ( = σ 2 I) and movement of the optimum along a single dimension (v = (v 1 , 0, . . . , 0) ). In this transformed phenotype space, the effects of selectional and mutational correlations are entirely captured by the M-matrix (and, in particular, the orientation of its leading eigenvector/first principal component) relative to the direction of environmental change. Furthermore, all vectors (e.g., z, v, α) are measured relative to the average width of the distribution of new mutationsm.

THE ADAPTIVE-WALK APPROXIMATION
The aim of this article is to investigate the distribution of adaptive substitutions φ(α), that is, the distribution of the effects of those mutations that eventually go to fixation and contribute to adaptation. Our main analytical tool will be the "adaptive-walk approximation." Following Kopp and Hermisson (2009b), this approximation is based on the simplifying assumption that whether a new beneficial mutation goes to fixation or is lost by drift is determined immediately after its appearance and that, in the former case, fixation occurs instantaneously. Therefore, the population can be considered monomorphic nearly all of the time, and adaptation occurs as a series of discrete "steps," which together will be referred to as an "adaptive walk" (Kauffman 1993;Orr 2000). This approximation ignores interactions between cosegregating mutations, such as epistasis, linkage, and Hill-Robertson interference (Hill and Robertson 1966).
Adaptive walks can easily be simulated using the following algorithm: (1) draw the waiting time for a new mutation from an exponential distribution with parameter /2 (where is a standard measure for the population-and genome-wide mutation rate); (2) draw the size of the mutation from its distribution p(α) (eq. 4); (3) accept the mutation (i.e., perform an adaptive step) with its fixation probability p fix (x, y, t) ≈ 2s (x, y, t) for s (x, y, t) > 0 0 f o rs (x, y, t) ≤ 0 (6) (Haldane 1927), where y is the current population phenotype, x = y + α is the mutant phenotype, and denotes the selection coefficient of the mutant x in a wild-type population with phenotype y at time t. Note that equation (6) assumes s to be small and neglects chance fixations of deleterious mutations. In some simulations, we also used the slightly more accurate approximation p fix ≈ max(0, 1 − exp(−2s)). (Even more accurate approximations exist that account for the change in the selection coefficient during the fixation process due to the environmental change, see Uecker and Hermisson 2011; however, within the simple framework of the adaptive-walk model we do not obtain further improvement.) The distribution of adaptive substitutions In the following we derive an analytical expression for the distribution φ(α|y) of the size α = x − y of the next adaptive substitution given an initial phenotype y at time t = 0. First, equation (7) can be approximated by with (provided λ x,y = 0). That is, with a linearly moving optimum, the selection coefficient increases or decreases approximately linearly over time, where λ x,y is the rate of change and τ x,y is the time when s reaches zero (the "lag time" in the terminology of Kopp and Hermisson 2007). This time dependence of the selection coefficient is illustrated in Supporting Information 1. The distribution φ(α|y) can then be calculated in four steps.
The instantaneous rate of substitutions. We denote by g(t, y) the rate at which substitutions of any kind happen at time t. g(t, y) is given by where χ(t, y) = {x | s (x, y, t) > 0} is the set of all mutant phenotypes with positive selection coefficient at time t. The integrand in (9) is simply the product of the probability that a mutation with phenotype x arises ( p(α)/2) and its probability of fixation approximated as (2s(x, y, t), see eq. 6).

The waiting-time distribution
We denote by F(t|y) the probability that no fixation has happened before time t (thus, 1 − F(t|y) is the cumulative distribution function of the waiting time to the next fixation). From the theory of Poisson processes, The conditional distribution of step sizes The distribution of step sizes, given that the step occurs at time t, is simply proportional to the distribution of new mutations weighted by the selection coefficient (Gillespie 1983;Kopp and Hermisson 2009b): The distribution of step sizes Finally, the unconditional distribution of the size of the next adaptive step can be calculated by integrating over all possible waiting times (see Kopp and Hermisson 2009b), yielding where f (t|y) = (1 − F(t|y)) is the density of the waiting-time distribution. Equation (12a) can also be written as where τ x,y is given by equation (8c).
The parameter γ Supporting Information 3 shows that the distribution of step sizes in the adaptive-walk approximation depends only on the distribution of new mutations and the composite parameter where the term in the numerator can be interpreted as the rate of environmental change relative to the width of the fitness landscape in the direction of the moving optimum. If selection is isotropic, equation (13) reduces to which is equivalent to the γ defined by Kopp and Hermisson (2009b) for the single-trait case, except for differences in notation, and is independent of n. Here,σ/m describes the mean width of the fitness landscape relative to the mean effect size of new mutations, and (σ/m) −2 can be seen as a scale-free measure for the strength of stabilizing selection. The product of this term and the population-wide mutation rate determines the "adaptive potential" of the population, γ can, thus, be interpreted as a scaled rate of environmental change (how fast the population needs to adapt relative to how readily it can adapt). In particular, it can be used to distinguish two limiting cases (Kopp and Hermisson 2009b).
If γ is small, the population can easily follow the optimum. The adaptive process is, therefore, environmentally limited, and the distribution of adaptive substitutions is primarily determined by the lag time τ x,y , which determines when a mutation of effect α = x − y becomes beneficial ("dynamic sieve" sensu Kopp and Hermisson 2009b). If γ is large, the population will follow the optimum with a large lag. In this case, the adaptive process is genetically limited, and the distribution of adaptive substitutions is largely determined by the distribution of new mutations p(α) ("static sieve" sensu Kopp and Hermisson 2009b). Numerical values for γ in these two regimes are discussed in Supporting Information 3.

The environmentally limited regime
In the environmentally limited regime, the Gaussian distribution of new mutations, p(α), can be approximated by a uniform distribution p u (α) with the same density at α = 0, that is, (see Kopp and Hermisson 2009b). This approximation is justified if the optimum moves so slowly that all beneficial mutations are small (α close to 0). It allows to directly calculate several properties of the distribution of adaptive substitutions. In particular, if the wild-type phenotype y = 0 (i.e., the population is perfectly adapted at time t = 0), the distribution of the "first" substitution (and all its moments) can be calculated analytically (Supporting Information 4).

INDIVIDUAL-BASED SIMULATIONS
In addition to our adaptive-walk approximation, we conducted individual-based simulations (implemented in C++, available at doi:10.5061/dryad.534f0; see Bürger 2000; Kopp and Hermisson 2009b), which allow multiple mutations to segregate simultaneously, while making additional assumptions about the genetic architecture of the selected traits, the life cycle of individuals and the regulation of population size. The simulations follow the evolution of a population of individuals with discrete and nonoverlapping generations. Individuals are characterized by L unlinked diploid loci, which additively determine the n-dimensional phenotype z. According to the universal-pleiotropy assumption, each allele at each locus is specified by a vector of contributions to the n traits. We neglect environmental variance and, therefore, equate genotypic and phenotypic values. Mutations occur at rate u per (diploid) locus and have effects drawn from the distribution p(α) (eq. 4). Each generation, the following steps are performed: (1) Viability selection: Individuals are removed with probability 1 − w(z) (eq. 1).
(2) Population regulation: If, after selection, the population size N exceeds the carrying capacity K , N − K randomly chosen individuals are killed (Bürger 2000).
(3) Reproduction: The surviving individuals are randomly assigned to mating pairs, and each mating pair produces exactly B offspring (typically, B = 4). Note that, with this procedure, the effective size of a well-adapted population exceeds the census size (e.g., for B = 4, N e = 4/3N ; Bürger 2000, p. 274). The offspring genotypes are derived from the parent genotypes by taking into account segregation, recombination, and mutation.
To monitor adaptive substitutions, the program keeps track of the genealogical relationship between the alleles at a given locus. A substitution is recorded whenever there is a change in the root of such an "allele tree" (i.e., when the surviving alleles get a new most recent common ancestor). This is equivalent to calling an allele fixed if the entire population has been taken over by that allele or its descendants (e.g., Gillespie 1993;Park and Krug 2007).
In all simulations, the initial population contained K = 1000 identical, homozygous individuals with phenotype 0 (i.e., the population was perfectly adapted at time 0). The number of loci was set to L = 10 and the mutation rate per diploid locus to μ = 5 × 10 −6 per generation. This yields a population-and genome-wide mutation rate = 2N Lu = 0.2. We chose this value to limit complications from interference between alleles cosegregating at the same locus, which have been thoroughly studied for the one-trait case in Kopp and Hermisson (2009b) (see Discussion). When comparing individual-based simulations to adaptive-walk simulations with differing , the speed of environmental change v was adjusted accordingly to reach the same value of γ (eq. 13). Simulations were stopped after 1000 substitutions had been recorded. Alternatively, we only recorded the first adaptive substitution for 1000 replicate runs to study the initial phase of the adaptive process. Finally, for some parameter combinations, the simulations terminated because the population went extinct (e.g., if the environmental change was too fast).

Results
Our primary interest is the distribution of adaptive substitutions, that is, the distribution φ(α) of the effects of those mutations that go to fixation and contribute to adaptation (eq. 12).

AND PHENOTYPIC COMPLEXITY
Key properties of the distribution of adaptive substitutions can already be seen from a simplified model in which both mutation and selection are isotropic. Note that any model in which the two matrices M and are proportional to each other, that is, have the same shape and orientation, can be reduced to this case via the transformation described in Supporting Information 2. The same holds true for any model under the environmentally limited regime, in which the shape of the M-matrix is irrelevant (because any distribution of new mutations can be approximated by a uniform distribution).
In the isotropic model, the distribution of adaptive substitutions is symmetric around the direction of the moving optimum. Figure 1 shows this distribution in adaptive-walk simulations with n = 2 traits. The marginal distribution in the direction of the optimum has an intermediate mode and resembles a gamma distribution, in accordance with previous results for the one-dimensional moving-optimum model (Kopp and Hermisson 2009b), and in contrast to the exponential pattern predicted for the classical Fisher model with constant selection (Orr 1998). Although the population always follows the optimum, pleiotropic side effects of fixed mutations frequently lead to maladaptation of the traits under pure stabilizing selection. The distribution of these deviations is bell-shaped and centered around zero (Fig. 1).
For small γ, explicit analytical results can be obtained for the distribution of the first adaptive substitution in the environmentally limited regime (Supporting Information 4, eq. S19). In particular, assuming v = (v 1 , 0, . . . , 0), the mean and variance in the direction of the optimum are given by where η(n) = π n 2 (n+3) (2+ n 2 ) and (•) denotes the gamma function. Interestingly, the coefficient of variation √ Var(α 1 |0)/E(α 1 |0) depends only on n (see Fig. S4 2). The variance in directions orthogonal to the optimum is given by Additional results regarding higher moments of α 1 and α 2 , the total step size α and the magnitude of pleiotropic deviations are given in Supporting Information 4.
In accordance with previous findings (Collins et al. 2007;Hermisson 2007, 2009a, b), equations (16)-(18) show that the mean step size in the direction of the optimum increases with the scaled rate of environmental change γ, and so does the magnitude of pleiotropic deviations (Figs. S4 1, S5 1). This fundamental relationship also holds true over the entire adaptive walk and beyond the environmental limit (Figs. 2, S5 2, S5 3).
Some discussion is warranted regarding the dependence of the mean step size on the average variance of mutational effects m 2 . Increasingm 2 decreases γ and, consequently, leads to a reduced mean step size in the transformed phenotype space (see Supporting Information 2, eq. S28), where phenotypes are measured relative tom. When phenotypes are measured in arbitrary units, however, the mean step size increases withm (eq. 16, see also Kopp and Hermisson 2009b). The reason is that an increase inm reduces the rate of appearance of small mutations (and, hence the parameter p 0 in the environmental limit, see eq. 15), which reduces the ability of the population to follow the optimum closely.
A key result of our analysis is that, for a given speed of environmental change, the mean step size in direction of the optimum increases with the number of traits under selection, that is, with the level of pleiotropy or organismic complexity (eq. 16, Figs. 2, S4 1, see also Figs. S5 2 and S5 3), and a similar result also holds for fitness (Fig. S5 4). At first, this result seems to contradict the "cost of complexity" argument from Fisher's model, which states that, in complex organisms, large mutations are unlikely to contribute to adaptation. The explanation is that, precisely because fewer mutations are beneficial if n is large (because there are more directions in which they can "go wrong"; Orr 1998), the time to the first step increases (see Supporting Information 4). By this time, the optimum has already moved considerably, such that also large mutations are beneficial (see Supporting Information 1 and Fig. S1 1), even if they have significant pleiotropic effects. These effects-in particular, the increased waiting time between adaptive substitutions-also affect population persistence: as shown in Figure S5 5, the mean time to extinction decreases with organismic complexity, and so does the maximal rate of environmental change a population can tolerate (Bürger and Lynch 1995).

SELECTIONAL AND MUTATIONAL CORRELATIONS
To study the orientation of the distribution of step sizes in the n-dimensional phenotype space, we now consider a model with correlated selection and correlated mutations. We will assume that the angle between the direction of the optimum v and the leading eigenvector of and/or M is 45 • . More precisely, the optimum moves along the first trait axis (v = (v 1 , 0, . . .) ), whereas the fitness landscape and/or the distribution of new mutations are concentrated along the main diagonal: All diagonal elements of (M) are equal to σ 2 (m 2 ) and all off-diagonal elements have magnitude ρ σ 2 (ρ M m 2 ), where 1 > ρ ≥ 0 (1 > ρ M ≥ 0 ) is the magnitude of selectional (mutational) correlation. In this case, the fitness landscape (distribution of new mutations) is symmetric around the leading eigenvector of (M), ν = (1, 1, 1, . . .). We first study the effects of mutational and selectional correlation separately. Exemplary adaptive walks for strong correlations are shown in Figure 3. Figure 4 shows the multivariate distribution of adaptive substitutions, φ(α), for different strengths of selectional and mutational correlations under varying speeds of environmental change for n = 2 traits. As in the isotropic case (Fig. 1), the distribution φ(α) is biased toward the direction of the optimum, with pleiotropic side effects of fixed mutations on average being neutral (Figs. S5 6, S5 7). The shape of the distribution, however, critically depends on the interaction between the type and strength of correlations and the rate of environmental change. Mutational correlations tend to align the distribution of adaptive substitutions along the leading eigenvector of M, with stronger mutational correlations leading to stronger correlation in step sizes (Figs. 4,5,and S5 8 top left). This effect is strongest in fast-changing environments and gradually gets weaker as the rate of environmental change decreases (Fig. 5), until becoming almost unnoticeable. Selectional correlations similarly orientate the distribution of adaptive substitutions along to the leading eigenvector of (Figs. 4, 5 bottom left, S5 8). In contrast to mutational correlations, however, their impact is strongest if environmental change is slow (for small γ and the first step, the correlation is given by ρ Var (α2|0) Var(α1|0) ≈ ρ , see eq. S36). Correlations in step sizes remain almost unchanged for a broad range of rates v 1 , before dropping off sharply once environmental change gets sufficiently fast.
These results still hold true when mutational and selectional correlations are both present but with opposite signs. As shown in Figure 6, the correlations in step sizes resemble the selectional correlations if environmental change is slow and resemble the mutational correlations if environmental change is fast. At intermediate rates of environmental change, the two effects cancel, and correlations in step sizes are close to zero.
Mutational and selectional correlations depend on the coordinate system in which multivariate phenotypes are measured (i.e., on the definition of traits). As shown in Supporting Information 2, there is always a transformation to coordinates in which selection (but not necessarily mutation) is isotropic. The key question, therefore, is whether the distribution of new mutations is aligned with the fitness landscape (in terms of the eigensystems/principle components of the matrices M and ). Our results can, thus, be reformulated as follows: The distribution of adaptive substitutions reflects the shape and orientation of the fitness landscape if adaptation is environmentally limited (i.e., if the optimum changes slowly), whereas it mirrors the distribution of new mutations (but with a mean shifted in the direction of the optimum) if adaptation is genetically limited (i.e., if the environment changes fast). Intuitively, as long as environmental change is slow, the population is close to the optimum and the shape of the distribution of new mutations is practically irrelevant, because only a small subset of new mutations from the center of their distribution can pass the selective sieve (Kopp and Hermisson 2009b). In contrast, if environmental change is fast, the population is far from the optimum, and the selective sieve has less impact on the adaptive process than the supply with new mutations. In the limit, pleiotropic side effects become negligible and the selection coefficient of new mutations depends only on their effect in the direction of the optimum.
Finally, mutational and selectional correlations also impact the trajectory of the mean phenotype (Figs. 3, 5, S5 8; see also Jones et al. 2004). In particular, strong mutational correlations can cause the mean phenotype to trail above and behind the moving optimum-an effect that has been phrased the "flying-kite effect" (Jones et al. 2004). Conversely, with strong (positive) selectional correlation, the phenotypic mean follows the optimum behind and below. In analogy to the flying-kite effect, we call this phenomenon the "diving-kite effect." Both effects can be explained by a deterministic model in which the change in the mean phenotype depends primarily on the leading eigenvector of M and the selection gradient β(t) (Fig. 7). Under strong mutational correlation, the change in mean phenotype is initially dominated by the leading eigenvector of M, causing the "rise of the kite," until it is balanced by the selection gradient pointing toward the optimum. Under strong selectional correlation, however, the selection gradient is initially orthogonal to the leading eigenvector of (i.e., the "ridge" of the fitness landscape), causing the mean phenotype to "dive." Again, the trajectory will gradually change until it is aligned with the direction of the moving optimum (where it is aligned with the axis of largest width of the fitness landscape). Observing either the flying or the diving kite requires a sufficiently fast-changing environment (the kite needs to be pulled strongly enough) and at least intermediate levels of mutational or selectional correlations (right column Figs. 5, S5_8). As the number of traits increases, the strength of both effects decreases on a per-trait basis, but their total strength increases (Supporting Information 4 and Fig. S4_1). Independently of the number of traits, the population on average takes smaller steps in the direction of the optimum as correlations (either selectional or mutational) become stronger (Figs. S5_2, S5_3).

ACCURACY OF THE APPROXIMATIONS
Our main analytical tool has been the adaptive-walk approximation with normally distributed mutational effects. When compared to individual-based simulations of an explicit genetic model, its performance is often surprisingly good (e.g., Figs. 4, 5, S5_2, S5_3, S5_8). Significant deviations occur mainly if the population-wide mutation rate is high ( 1), which, in violation of the adaptive-walk assumption, increases the probability of cosegregating beneficial mutations (Discussion and Figs. S5_9, S5_10). Because the adaptive-walk approximation does not account for population dynamics, it cannot be used to predict population persistence or extinction. Individual-based simulations show For slow environmental change, the normal distribution of new mutations can, furthermore, be approximated by an appropriate uniform distribution. The resulting approximation for adaptive walks works well for a broad range of small to intermediate rates of environmental change (see insets in Fig. S5_2). Naturally, this approximation cannot capture mutational correlations (see the poor fit for high values of ρ M and v 1 in Fig. S5_3). Note, however, that for sufficiently small rates of environmental change, mutational correlations can, indeed, be ignored (see above,Figs. 4,5).
Finally, we have attempted to approximate the distribution of adaptive substitutions over an entire adaptive walk by the distribution of the first step. This approximation works well in the one-trait case (Kopp and Hermisson 2009b), and in combination with a uniform distribution of new mutations, it is the only approach that allowed significant analytical progress (Supporting Information 4). With multiple traits, however, the first step makes a larger progress toward the optimum than the subsequent steps (Figs. 2, S5_2, S5_3). The reason is that the first step will always introduce maladaptive pleiotropic side effects, which become compensated for by subsequent substitutions. Some of these compensatory substitutions are "backward steps," which are beneficial, despite their effect being opposite to the direction of the optimum (e.g., Supporting Information 1, gray ellipse in Fig.  S1_1 and "backward steps" in Figs. 1 and 4). Consequently, the first-step approximation works less well as the number of traits increases (Fig. S5_1). Furthermore, with mutational or selectional correlations, the direction of the first step systematically deviates from the distribution of step sizes over the entire adaptive walk (see flying-and diving-kite effects above; for the case of selectional correlation, see Fig. S5_11).

Discussion
Environmental change forces populations to either adapt to the altered conditions or go extinct. In the absence of standing genetic variance for fitness, the outcome crucially depends on mutations, which provide the "genetic fuel" for adaptation, and selection that converts this resource into adaptive substitutions. Here, we have used analytical approximations and individual-based simulations to study the effects of pleiotropy or "organismic complexity" on the genetic basis of adaptation in gradually changing environments. In particular, we have investigated the distribution of adaptive substitutions (i.e., the distribution of the phenotypic effect sizes of fixed mutations) in populations following a moving optimum. Our results confirm and extend previous analysis of "adaptive walks" for single traits (Collins et al. 2007;Hermisson 2007, 2009a,b). We show that the distribution of adaptive substitutions is largely determined by a single composite parameter γ, which scales the rate of environmental change relative to the "adaptive potential" of the population and defines a continuum between environmentally and genetically limited adaptation (Kopp and Hermisson 2009b). In the environmentally limited regime (slow environmental change), the population follows the optimum closely, adaptive steps are small and their multivariate distribution mirrors the shape of the fitness landscape.
In the genetically limited regime, in contrast, the population follows the optimum with a large gap, adaptive steps are large and their distribution is determined primarily by the distribution of new mutations. We furthermore show that the mean effect size of fixed mutations increases with the degree of pleiotropy, in contrast to classical predictions from FGM under sudden environmental change. We now discuss these results in greater detail.

GENETICS OF ADAPTATION
In complex organisms, pleiotropy is widespread-that is, most mutations affect multiple traits simultaneously. Different traits, therefore, do not evolve independently (Lande 1979;Agrawal and Stinchcombe 2009;Walsh and Blows 2009). With this basic fact in mind, Fisher (1930) used his classical geometric model to argue for a predominance of small mutations in adaptive evolution. Although theoretical studies later pointed out that large beneficial mutations may, nevertheless, play an important role (Kimura 1983;Gillespie 1993;Orr 1998Orr , 2005a, they also confirmed that organisms pay a "cost of complexity" (Orr 1998(Orr , 2000Welch and Waxman 2003) in the form of a reduced rate of adaptation. With regard to individual substitutions, Orr (2000) found that more complex organisms make smaller steps when adapting toward a fixed optimum (with step size measured as the decrease in the absolute distance to the optimum, which is closely related to the fitness effect of a fixed mutation). This is in direct contrast to our results for a moving optimum, where increased complexity leads to larger step sizes, with respect to both phenotype and fitness (Figs. 2, S5_4). The main reason for this finding arises from the ecological differences between the classical FGM and the moving-optimum model (consequences of different mutation models are discussed below). In the classical Fisher model, the proportion of beneficial mutations decreases with organismic complexity. Thus, the more phenotypic traits, the longer one has to wait for a beneficial mutation to appear (as adding another trait adds yet another dimension where mutations can go wrong; eq. S18). Of course, this argument still holds true under a moving optimum. As more complex organisms have to wait longer for a beneficial mutation to appear, the optimum has already traveled farther, enabling larger mutations to become fixed. Thus, the moving-optimum model does not contradict the "cost of complexity" argument, but reveals yet another aspect of it.

SELECTIONAL CORRELATIONS
The impact of mutational and selectional correlations on the distribution of adaptive substitutions is a direct consequence of the general principle that the shape of this distribution depends on the scaled rate of environmental change (see above). In particular, if the rate of environmental change is slow, only mutations from the very center of the mutational distribution can pass the selective sieve (Kopp and Hermisson 2009b), making mutational correlations irrelevant relative to the shape of the fitness landscape. Conversely, if adaptation is genetically limited, the selective sieve has less impact on the adaptive process than the supply with new mutations. Between these two extremes, the distribution of adaptive substitutions will progressively take the orientation of the mutational distribution as the rate of environmental change increases (Fig. 6).
Our results reveal strong parallels between the distribution of adaptive substitutions and the evolution of the G-matrix describing standing genetic variation (see also below). Recent quantitative-genetic studies have shown that both mutational and selectional correlations can induce correlations in the G-matrix, under both constant stabilizing and moving-optimum selection (Jones et al. 2003(Jones et al. , 2004. A link between selectional correlation and genetic correlation has also been confirmed empirically (see Roff and Fairbairn 2012 for a recent meta-analysis). As shown in Figure S5_13, the distribution of adaptive substitutions closely matches the shape and orientation of the G-matrix. While this seems intuitive, it had not been shown by any previous study, and little is known about the relation between alleles in the standing variation and those that ultimately reach fixation (but see Hill 1982;Hill and Rasbash 1986a,b). This close correspondence between standing variation and fixed mutations might explain why the adaptive-walk approximation works surprisingly well even in populations with a high mutation rate (see Kopp and Hermisson 2009b). Quantitative-genetic studies have, so far, not systematically investigated how correlations in the G-matrix are affected by the rate of environmental change. It would be interesting to know whether the effects of mutational and selectional correlations on the G-matrix are similar to those on the distribution of adaptive substitutions.
Confirming previous results by Jones et al. (2004), our simulations showed that mutational and selectional correlations can cause systematic maladaption in traits under purely stabilizing selection (i.e., in directions orthogonal to the direction of the optimum). These "flying-" and "diving-kite" effects require that correlations are strong and the environmental change is sufficiently fast (i.e., in the genetically limited regime, see Figs. 3 and 5). Strong effects are, therefore, likely to be restricted to a narrow parameter range, where populations might often be on the brink of extinction.

FUTURE DIRECTIONS
Like all models, our study is based on a number of simplifying assumptions, which might constrain the generality of our results. In the following, we discuss the likely consequences of these assumptions, potential extensions of the model, and ways to test our predictions empirically. First, our model is based on the assumption of universal pleiotropy (Kacser and Burns 1981; for a review see Paaby and Rockman 2013). This assumption has been challenged recently, both because empirical levels of pleiotropy are rather low (median 1-7; Wang et al. 2010) and because true universal pleiotropy would induce unsustainably high costs ("the cost of complexity [...] should be more properly called the cost of pleiotropy"; Wagner and Zhang 2012, but see Hill and Zhang 2012a). Alternative approaches have, therefore, suggested modularity (Wagner and Altenberg 1996;Welch and Waxman 2003) or partial pleiotropy (Chevin et al. 2010;Lourenço et al. 2011) as a solution to this problem. Indeed, our model might be best interpreted as applying to a given module with a moderate level of pleiotropy. In any case, the relatively low number of traits assumed in most parts of this article is consistent with the degree of pleiotropy observed in natural populations (Martin and Lenormand 2006a,b;Wang et al. 2010). Thus, we expect our results to apply across a wide range of species facing environmental change.
Second, we assume the so-called Euclidean superposition model (Turelli 1985;Wagner 1988), where the distribution of mutational effects on a given trait is independent of complexity (see also Welch and Waxman 2003;Lourenço et al. 2011;Zhang 2012). Other studies (Orr 1998(Orr , 2000Wingreen et al. 2003; and an alternative model in Welch and Waxman 2003) have instead used a "constant total-effects model," in which the total mutational effect size ( α ) is constant across levels of complexity and, in consequence, the mean effect size on individual traits decreases. Indeed, this assumption explains part of our differences to Orr (1998Orr ( , 2000. More generally, it raises the question of which factor is more important in shaping the distribution of adaptive substitutions at different levels of complexity: the "pleiotropic scaling" of mutations (Wagner et al. 2008) or the mode of environmental change. To address this issue, we conducted additional simulations, which combined constant and moving-optimum selection with the Euclidean superposition and constant total-effects models. These simulations yielded three main results (Fig. S5_12). First, the moving-optimum model behaves qualitatively similarly under both mutation models; in particular, average step size increases with complexity (at least as long as adaptation remains environmentally limited). This shows that our main results are robust to considerable variation in the pleiotropic scaling of mutational effects. Second, with a constant optimum, the mutation model does make a qualitative difference for the mean step size in the direction of the optimum (but not for total step size α ), which decreases with complexity under the constant total-effects model, but increases under the Euclidean superposition model. Third, a fundamental difference between constant and movingoptimum selection, which is independent of the mutation model, is seen at the level of the selection coefficients of fixed mutations, which decrease with complexity under a constant optimum but increase with complexity under a moving optimum. In summary, the mode of environmental change plays a fundamental role in shaping the distribution of adaptive substitutions and, in many cases, overrides the effects of pleiotropic scaling. Nevertheless, a better understanding of pleiotropic scaling-both empirically (Hermisson and McGregor 2008;Wagner et al. 2008;Wang et al. 2010;Wagner and Zhang 2011;Hill and Zhang 2012b) and with respects to its theoretical consequences-clearly is an important topic for future research.
Third, our adaptive-walk approximation assumes that evolution proceeds as a series of (hard) selective sweeps originating from new mutations. This follows the tradition of models based on a strong-selection-weak-mutation approximation (Gillespie 1983;Orr 1998Orr , 2005a. In contrast, quantitative-genetic models assume that virtually all adaptation stems from standing genetic variation (in the context of the moving-optimum model, see, e.g., Bürger and Lynch 1995;Gomulkiewicz and Holt 1995;Jones et al. 2004Jones et al. , 2012Zhang 2012;Chevin 2013;see also Fig. S5_13), and the importance of standing variation is well documented empirically (Hermisson and Pennings 2005;Barrett and Schluter 2008;Gomulkiewicz and Houle 2009;Teotónio et al. 2009;Jerome et al. 2011;Domingues et al. 2012;Messer and Petrov 2013). So far, quantitative-genetic models with a multidimensional moving optimum have focused either on the risk of population extinction (Gomulkiewicz and Houle 2009) or the maintenance and structure of genetic variation (Jones et al. 2004(Jones et al. , 2012. To our knowledge, very little is known about the distribution of phenotypic effect sizes of adaptive substitutions when adaptation occurs from standing genetic variation. Because in this case the adaptive process is likely to have very different properties (e.g., adaptation could be faster with on average smaller mutations becoming fixed Barrett and Schluter 2008;Rockman 2012), this should be an important topic for future research.
Fourth, in accordance with the adaptive-walk approximation, most of our simulations (including those in Fig. 6) assumed a relatively (but not unrealistically) low population-wide mutation rate . As such, we ignore effects of interactions between cosegregating beneficial mutations. For the one-dimensional case, Kopp and Hermisson (2009b) showed that high mutation rates in combination with low recombination (or a small number of loci) lead to an increase in the mean size of adaptive substitutions, due to Hill-Robertson interference (Hill and Robertson 1966;Gerrish and Lenski 1998). At high recombination rates (or with a large number of unlinked loci), in contrast, the mean step size decreases as a result of epistasis for fitness (due to stabilizing selection). Individual-based simulations suggest that these results also hold true for the multivariate case (Figs. S5_9, S5_10). Note, however, that the strength of interference is expected to decrease with increasing complexity, because the rate of beneficial mutations decreases.
In addition, Hill-Robertson interference also influences how the distribution of adaptive substitutions is affected by mutational and selection correlations. In particular, correlations between adaptive substitutions increase with increased linkage in the presence of mutational correlations (Fig. S5_9), but decrease with linkage in the presence of selectional correlations (Fig. S5_10). Thus, at high mutation rates, increasing linkage has a similar effect as increasing the scaled rate of environmental change γ. This makes intuitive sense, because interference weakens the efficiency of selection (Gerrish and Lenski 1998;Weissman and Barton 2012), which brings the adaptive process closer to the genetically limited regime.
Finally, our adaptive-walk approximation does not consider population dynamics. By setting an arbitrary extinction threshold with respect to mean fitness, we found that the maximal rate of environmental change a population can tolerate (Bürger and Lynch 1995), as well as the mean time to extinction, decreases with the number of traits ( Fig. S5_5; this result is also supported by individual-based simulations). In particular, in complex organisms, long-term persistence in the face of an indefinitely moving optimum seems to be possible only in the environmentally limited regime (γ 0.1), that is, when adaptation is not limited by the availability of new mutations. In the genetically limited regime, in contrast, populations can only persist for a limited amount of time (e.g., fast environmental change followed by a period of stasis). Over shorter time scales and with high mutation rates (large , as in Jones et al. 2004), population persistence can also be facilitated by adaptation from standing genetic variation (Bürger and Lynch 1995;Barrett and Schluter 2008;Gomulkiewicz and Houle 2009). Indeed, our parameter γ is structurally similar to expressions describing the equilibrium phenotypic lag in quantitative-genetic models of adaptation to a moving optimum (eq. 5 in Jones et al. 2004; see also eq. 8a in Bürger and Lynch 1995). The interpretation is analogous: the lag increases with the speed of environmental change, and decreases with the strength of selection and the amount of (standing) genetic variation.
Our model makes a number of concrete predictions ( Table 1) that can be tested empirically, even though such tests will certainly be challenging. The most direct approach is experimental evolution (for reviews see Elena and Lenski 2003;Kawecki et al. 2012;Barrick and Lenski 2013). Although the majority of studies have employed constant conditions (Reusch and Boyd 2013), Collins (2011b) recently urged for more studies in gradually changing environments. Microorganisms such as bacteria, yeast, or algae can be cultivated in media where an environmental factor such as temperature (Hietpas et al. 2013), salinity (Bell and Gonzalez 2009;Lachapelle and Bell 2012), pH (Hughes et al. 2007), the availability of nutrients (Collins 2011a), or the concentration of stressors such as antibiotics (Perron et al. 2008;Lindsey et al. 2013) or pollutants (Adamo et al. 2012) is gradually changed. Until now, these studies were mainly used to investigate the probability of "evolutionary rescue" (Gonzalez et al. 2013). Recent advances in sequencing technologies (reviewed in Metzker 2010), however, make it possible to conduct real-time genome-wide analyses and to map genetic changes to their effects on phenotype and fitness (Barrick et al. 2009;Barrick and Lenski 2013), such that the distribution of adaptive substitution over entire adaptive walks can be analyzed.
In natural populations, where only present-day data are available, the most promising approach for studying the genetic basis of adaptation is the analysis of quantitative-trait loci (QTLs) in diverging populations. For example, Langlade et al. (2005) identified QTLs for leaf shape in two species of Antirrhinum and postulated a sequence of substitutions that can traverse the "allometric space" between them. Albert et al. (2008) and Rogers et al. (2012) analyzed genetic differences between ancestral marine and derived freshwater populations of sticklebacks (Gasterosteus aculeatus). Rogers et al. (2012) compared two sets of freshwater populations and showed that those whose environment (and presumably phenotypic optimum) is more different from that of the marine populations (with respect to salinity and presence of predators) displayed a higher frequency of large-effect QTLs. Albert et al. (2008) crossed an ancestral Pacific and a highly derived benthic freshwater form and found a gamma-like The distribution of adaptive substitutions is approximately exponential with respect to phenotype and fitness. Accordingly, most fixed mutations are of small effect and only a few large-effect alleles become fixed when approaching the constant optimum. The farther the optimum is away (i.e., the harsher the sudden environmental change) the larger the mutational effects that get fixed.

Gradual change
The distribution of adaptive substitutions with respect to phenotype (in the direction of the optimum or total effect) and fitness is gamma-like with an intermediate mode (Figs. 1, S5_4). Thus, when following the moving optimum, most adaptive substitutions are of intermediate effect with only a few large-effect alleles becoming fixed. Scaled rate of environmental change The faster the rate of environmental change relative to the adaptive potential, the larger the mutational effects that become fixed. Holds true with respect to phenotype and fitness (Figs. 2, S5_4). With increasing rate of environmental change the distribution of fitness effects becomes more asymmetric.

Complexity/Pleiotropy
Mean effect of adaptive substitutions with respect to phenotype and fitness increases as the number of traits affected by a single mutation increases (Figs. 2, S5_12).

Mutational correlation
If the rate of environmental change is fast, the distribution of adaptive substitution mirrors the mutational distribution (Fig. 6).

Selectional correlation
If the rate of environmental change is slow, the distribution of adaptive substitution reflects the shape of the fitness landscape (Fig. 6).
distribution of QTL effect sizes with an intermediate mode. Taking into account the detection limits for small-effect QTLs (Otto and Jones 2000), they interpreted this result as support for FGM with constant selection (i.e., the difficulty in identifying small QTLs would turn the predicted exponential distribution into an observed gamma-like distribution). However, both studies could, in principle, also be interpreted as showing the outcome of adaptation to a moving optimum (as briefly discussed in Schluter et al. 2010;Rogers et al. 2012), which directly predicts a distribution of effect sizes with an intermediate mode. More stringent tests of the present theory would require studying populations for which a moving optimum can be assumed a priori (e.g., comparisons of microalgae from pristine habitats with populations known to have experienced gradual eutrophication). Even then, the difficulty in detecting small-effect substitutions will remain a major challenge (Otto and Jones 2000).

CONCLUSION
Natural populations are constantly forced to adapt to changing environments, a process that takes place in a high-dimensional phenotype and genotype space. Along with previous studies, our analysis of the moving-optimum model shows that the genetic basis of this process depends critically on the tempo and mode of environmental change. In particular, our environmentally and genetically limited regimes lead to qualitative differences in the distribution of adaptive substitution, with respect to its mean, shape, and correlation patterns. Long-term persistence is likely restricted to the environmentally limited regime-where adaptation proceeds "smoothly" in small steps-but the parameter range for this regime is reduced in complex organisms.

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's website: Supporting Information 1. The selection coefficient. Figure S1_1. The time dependence of selection in the two-dimensional moving-optimum model. Supporting Information 2. Transformation of phenotype space. Supporting Information 3. The parameter γ. Supporting Information 4. The environmentally-limited regime. Figure S4_1. Dependence of various components of the first adaptive stepα in the transformed phenotype space, as a function of the scaled rate of environmental change γ for various numbers of traits n, assuming a uniform distribution of new mutations. Figure S4_2. Coefficient of variation of the size of the first step in the direction of the moving optimum, √ Var(α 1 )/E(α 1 ). Supporting Information 5. Supplementary figures. Figure S5_1. Distribution of the size of the first adaptive substitution for two (top row) and four (bottom row) traits and two different rates of environmental change v 1 . Figure S5_2. Mean sizeα of adaptive substitutions in the direction of the moving optimum for the first step (lines) and over the entire adaptive walk (symbols), as a function of the rate of environmental change v 1 , for n = 2 and three traits and different values of selectional correlation ρ . Figure S5_3. Mean sizeα of adaptive substitutions in the direction of the moving optimum, for the first step (lines) and over the entire adaptive walk (symbols), as a function of the rate of environmental change v 1 , for n = 2 and three traits and different values of mutational correlation ρ M . Figure S5_4. Mean selection coefficient of adaptive substitutions as a function of the rate of environmental change for various strengths of mutational (dashed lines) and selectional (solid lines) correlations ρ M and ρ . Figure S5_5. The mean time to extinction for different levels of phenotypic complexity n as a function of the rate of environmental change. Figure S5_6. The multivariate distribution of the first adaptive substitution (left) and over the entire adaptive walk (right) for n = 2 traits, when the optimum moves slowly in the direction of the first trait and the effects of new mutations are strongly correlated (ρ M = 0.9). Figure S5_7. The multivariate distribution of the first adaptive substitution (left) and over the entire adaptive walk (right) for n = 2 traits, when the optimum moves slowly in the direction of the first trait and selection is strongly correlated (ρ = 0.9). Figure S5_8. The impact of mutational and selectional correlations on the distribution of adaptive substitutions for n = 3 traits. For details, see Figure 5 of the main text. Figure S5_9. The effects of linkage and interference between cosegregating alleles on the mean step size in direction of the moving optimum α 1 (left) and the correlation between adaptive substitutions ρ α 1 ,α 2 (right) under strong mutational correlations ρ M = 0.9. Figure S5_10. The effects of linkage and interference between cosegregating alleles on the mean step size in direction of the moving optimum α 1 (left) and the correlation between adaptive substitutions ρ α 1 ,α 2 (right) under strong selectional correlations ρ = 0.9. Figure S5_11. The multivariate distribution of the first adaptive substitution in adaptive-walk simulations with strong selectional correlation (ρ = 0.9), illustrating the diving-kite effect (the negative bias in the α 2 -direction) present for fast (v 1 = 0.1) but not for slow (v 1 = 10 −5 ) environmental change. Figure S5_12. Comparison of the first adaptive substitution under the constant total effect (darker gray) and the Euclidean superposition model (lighter gray) for Fisher's geometric model with constant selection (left) and a moving optimum (right). Figure S5_13. The G-matrix (gray ellipse) and the 90% confidence ellipse of the distribution of adaptive substitutions (dark ellipse) under strong mutational (ρ M = 0.9) and selectional (ρ = 0.9) correlation for various rates of environmental change (v 1 ).