# Modeling of autophagy-related gene expression dynamics during long term fasting in European eel (Anguilla anguilla)

Valérie Bolliet, Jacques Labonne, Laure Olazcuaga, Stéphane Panserat, Iban Seiliez

► To cite this version:

Valérie Bolliet, Jacques Labonne, Laure Olazcuaga, Stéphane Panserat, Iban Seiliez. Modeling of autophagy-related gene expression dynamics during long term fasting in European eel (Anguilla anguilla). Scientific Reports, 2017, 7 (1), 10.1038/s41598-017-18164-6. hal-01695895

HAL Id: hal-01695895

https://hal.science/hal-01695895

Submitted on 29 Jan 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_2.png" /></figure>

OPEN

# Modeling of autophagy-related gene expression dynamics during long term fasting in European eel (Anguilla anguilla)

Received: 13 September 2017

Accepted: 7 December 2017

Published online: 20 December 2017

Valérie Bolliet<fnref n="1" />,², Jacques Labonne<fnref n="1" />,², Laure Olazcuaga<fnref n="1" />,², Stéphane Panserat³,⁴ & Iban Seiliez ID ³,⁴

Autophagy is an evolutionary conserved cellular self-degradation process considered as a major energy mobilizing system in eukaryotes. It has long been considered as a post-translationally regulated event, and the importance of transcriptional regulation of autophagy-related genes (atg) for somatic maintenance and homeostasis during long period of stress emerged only recently. In this regard, large changes in atg transcription have been documented in several species under diverse types of prolonged catabolic situations. However, the available data primarily concern atg mRNA levels at specific times and fail to capture the dynamic relationship between transcript production over time and integrated phenotypes. Here, we present the development of a statistical model describing the dynamics of expression of several atg and lysosomal genes in European glass eel (Anguilla anguilla) during long-term fasting at two temperatures (9 °C and 12 °C) and make use of this model to infer the effect of transcripts dynamics on an integrated phenotype – here weight loss. Our analysis shows long-term non-random fluctuating atg expression dynamics and reveals for the first time a significant contribution of atg transcripts production over time to weight loss. The proposed approach thus offers a new perspective on the long-term transcriptional control of autophagy and its physiological role.

Over the last decades, arrivals of European glass eels *Anguilla anguilla* dramatically dropped (by up to 90%), leading to an important decline in continental stocks of eels<fnref n="1" />. In the light of this decline, the European eel was included in 2007 in Appendix II of the Convention on International Trade of Endangered Species (CITES) and was listed in 2008 as critically endangered in the IUCN Red List of Threatened Species. A management framework for the recovery of the European eel stock was established in 2007 by the Council of the European Union. However, the restauration attempts are made difficult by gaps in knowledge of eel biology and the functioning of this panmictic population. Management of the population will benefit of a better understanding of the life cycle of eel and the underlying mechanisms.

The species life cycle spans over a large area and through various ecosystems from the Sargasso Sea, where spawning occurs, to the European continental shelf, where larvae undergo a first metamorphosis into glass eel². At this stage, glass eels display a large panel of migratory tactics, ranging from residency in marine water, to various degrees of upstream colonization through estuarine and freshwater ecosystems³. However, these different patterns of migration can have a strong impact on the fate of the population. Indeed, sex determination is environmental: individuals remaining downstream mostly develop in males and return earlier to the Sargasso Sea, whereas individuals colonizing upstream develop mainly into females and stay longer on the continent⁴,⁵. Therefore, variations in propensity to migrate upstream can have major consequences on the phenotypic structure (sex and size distribution) and thus on the persistence of the species. Such propensity to migrate has been proposed to be related to their energy status. Indeed, most glass eels starve during the estuarine migration⁶ and their ability to reach the upstream zones would closely depend on the mobilization of their energy stores. The

<footnote n="1">¹INRA, UMR 1224 ECOBIOP, F-64310 St Pée sur, Nivelle, France. ²Univ Pau & Pays Adour, UMR 1224 ECOBIOP, UFR Sciences et Techniques Côte Basque, Anglet, France. ³INRA, UMR 1419 Nutrition Metabolisme Aquaculture, F-64310 Saint Pée sur, Nivelle, France. ⁴Univ Pau & Pays Adour, UMR 1419 Nutrition Metabolisme Aquaculture, F-40000, Mont de Marsan, France. Valérie Bolliet and Jacques Labonne contributed equally to this work. Correspondence and requests for materials should be addressed to I.S. (email: iban.seiliez@inra.fr)</footnote>
propensity to migrate in glass eels may thus results from an adaptation of their metabolism and the establishment of a tightly controlled energy mobilizing system during fasting.

It is now widely accepted that macroautophagy (hereafter referred to as autophagy) plays a major role in mobilizing diverse cellular energy and nutrient stores, including proteins, carbohydrates and lipids during starvation in a large panel of taxa⁷. Autophagy is an evolutionarily conserved process in eukaryotes by which cytoplasmic cargo sequestered inside double-membrane vesicles is delivered to the lysosome for degradation. This 'self-eating' process not only rids the cell of intracellular misfolded proteins or damaged organelles but is also a key adaptive response to provide energy when nutrients are scarce. The level of autophagy increases during prolonged periods of cellular stress in many eukaryotic species (from yeast to mammals) in order to conserve energy and promote survival. In this context, and because glass eels undergo strong fasting while colonizing continental ecosystems, genetic- and environment-related variations of autophagy function can strongly affect their migratory ability as well as their sex and their subsequent survival. Therefore, characterization of functional variation of autophagy at an individual level is a key to address the evolutionary adaptive potential of this species in face of environmental changes.

Accumulating evidence indicates that the autophagic response to stress or starvation may proceed in two phases⁸. A rapid increase in the autophagic flux, which occurs within minutes or hours of exposure to stressful conditions and is entirely mediated by post-translational protein modifications, is generally followed by a delayed and protracted autophagic response that relies on the activation of specific transcriptional programs⁹. Such a biphasic mechanism is explained by the fact that under normal circumstances, most cells have relatively high amounts of autophagy-related proteins (ATG) readily available to generate autophagosomes for the first few hours of starvation¹⁰. But during prolonged starvation, with depletion of available ATG, the demand is accomplished by synthesizing new proteins through transcription and translation to maintain cellular homeostasis⁹.

The mechanisms underlying the regulation of autophagy have long attracted extensive concern of numerous scholars, and its acute transcriptional-independent regulation by nutrient-sensing signaling pathways is now well described. Thus, existing mechanistic knowledge has been used to build mathematical models that predict how the crosstalk between autophagy and other biological processes work together to determine cell fate decision in response to short-term stresses¹¹⁻¹⁵. Far less is known about the long-term transcriptional dependent regulation of autophagy. Indeed, although many signaling pathways and transcriptional factors involved in the transcriptional control of autophagy- and lysosome-related genes are now well described⁸, to our knowledge, no data are currently available on the dynamic aspects of this long-term process. Furthermore, while transcriptional regulation of autophagy genes has emerged as an important mechanism for ensuring the somatic maintenance and homeostasis during long period of stress¹⁶, the relationship between *atg* gene transcription and phenotypic changes remains largely unexplored both at cellular and tissue/organism levels. In our opinion, the key to make progress in this very field is to actually recognize that integrated phenotypes are complex objects that result from various processes which effects accumulate and interact over time. Therefore, our ambition is to target at long-term/temporal dynamics of transcripts to understand the construction of integrated phenotypes – here weight loss.

In the present study, we propose a statistical exploration of the expression of key autophagy- (*ulk1*, *lc3b*, *atg7*, *atg12*) and lysosome- (*cathepsin A*, *cathepsin D*, *cathepsin F*, *cathepsin L*) related genes in glass eel during long-term fasting (60 days). These genes were chosen based on their reported significant role in the biogenesis of autophagosomes¹⁷ and the protease activity of lysosomes¹⁸⁻²⁰, which is recognized to be responsible for the degradation of muscle proteins during fish migration, maturation and starvation²¹. We analyze the general temporal patterns of expression throughout time to detect linear and non-linear trends, and we evaluate inter-individual variation in gene expression. Then, we extend the analysis by testing the hypothesis that autophagy related transcripts produced over time can contribute to changes in an integrated phenotype (the weight loss), and by assessing how these genes interact with markers known to drive energy supply through catabolism of proteins and lipids (i.e., *glutamate dehydrogenase (gdh)* and *carnitine palmitoyltransferase I (cpt1)*, respectively).

In the coming decades, predicted change of rainfall patterns throughout seasons are expected to modify the thermal dynamics in rivers and estuaries²². Metabolic demand is increasing with temperature and we hypothesize that autophagy will also increase in glass eels to provide energy for migration. The ability to increase the expression of autophagy- and lysosome-related genes under increasing temperature may thus be critical for migration. In this context, two relatively contrasting temperatures (9 °C and 12 °C) belonging to a natural range are also challenged in the above presented analyses.

## Results

**Gene expression dynamics.** Since gene expression measurements require the sacrifice of glass eels, we only had access to punctual gene expression (at the body scale) at the time of death for each individual. Therefore, the statistical approach that follows aimed at using population scale (or group scale) data to infer possible individual dynamics.

We assumed that temporal dynamics of gene expression between individuals of the same species – and the same population by virtue of sampling – share to some extent some similarity when submitted to the same environmental stressors – here, fasting and temperature. Under such hypothesis, since glass eels were sampled at different fasting duration, any non-random pattern in the temporal dynamics of gene expression of the whole population can inform on the temporal dynamics of each individual profile. But non-conversely, it is paramount to realize that a random temporal pattern at the population scale can either indicate fully uncorrelated non-random profiles between individuals, or fully correlated random profiles between individuals, although the latter is far less likely.
<table><caption>**Table 1.** Gene expression models selection. For each gene, the Akaike Information Criterion (AIC) has been calculated; the minimal value is underlined, indicating the best trade-off between parsimony and fit. The null model represents the mean of data, the simplest possible model.</caption><thead><tr><th colspan="10">Genes</th></tr><tr><th>Models</th><th>atg7</th><th>atg12</th><th>lc3b</th><th>ulk1</th><th>catha</th><th>cathd</th><th>cathf</th><th>cathl</th><th>gdh</th><th>cpt1</th></tr></thead><tbody><tr><td>Null model</td><td>198.4091</td><td>200.9975</td><td>195.0637</td><td>274.8661</td><td>217.5909</td><td>303.4851</td><td>299.9847</td><td>288.2715</td><td>211.8984</td><td>332.8236</td></tr><tr><td>Temperature only model</td><td>200.3295</td><td>202.651</td><td>192.973</td><td>271.6297</td><td>219.0855</td><td>305.1524</td><td>301.9706</td><td>286.9918</td><td>212.971</td><td>325.722</td></tr><tr><td>Linear model</td><td>188.2157</td><td>168.6486</td><td>162.3492</td><td>272.4148</td><td>195.9116</td><td>297.8339</td><td>298.2102</td><td>258.8775</td><td>211.7331</td><td>232.7246</td></tr><tr><td>2nd order polynomial</td><td>188.0247</td><td>142.2302</td><td>147.0513</td><td>240.335</td><td>187.7946</td><td>297.0104</td><td>296.297</td><td><strong>243.8806</strong></td><td>212.6022</td><td>216.468</td></tr><tr><td>3rd order polynomial</td><td>189.7838</td><td>130.495</td><td>142.1297</td><td>239.5692</td><td>188.6795</td><td>295.394</td><td>296.4119</td><td>244.4375</td><td>202.7846</td><td>212.8415</td></tr><tr><td>4th order polynomial</td><td>191.3504</td><td>125.182</td><td>143.1294</td><td>239.5406</td><td>187.5967</td><td>297.3427</td><td>298.4112</td><td>245.4981</td><td>199.1094</td><td>213.5762</td></tr><tr><td>5th order polynomial</td><td>185.5638</td><td>126.3158</td><td>141.1464</td><td>241.3106</td><td>183.3651</td><td>298.9851</td><td>300.4111</td><td>245.94</td><td>199.8602</td><td>215.2168</td></tr><tr><td>6th order polynomial</td><td><strong>177.9332</strong></td><td>126.7975</td><td><strong>135.2623</strong></td><td>240.2511</td><td>176.1458</td><td>295.1183</td><td><strong>294.0369</strong></td><td>247.1009</td><td>200.1375</td><td>213.4855</td></tr><tr><td>Linear model + temperature</td><td>189.6107</td><td>163.6934</td><td>164.0462</td><td>271.2542</td><td>197.3702</td><td>297.8925</td><td>299.7345</td><td>260.7578</td><td>213.4121</td><td>234.1897</td></tr><tr><td>2nd order polynomial + temperature</td><td>189.713</td><td>139.9847</td><td>149.0366</td><td>241.0182</td><td>188.1298</td><td>297.6679</td><td>298.145</td><td>245.7898</td><td>214.0683</td><td>218.4512</td></tr><tr><td>3rd order polynomial + temperature</td><td>191.5685</td><td>130.6108</td><td>143.9028</td><td>240.8709</td><td>188.2381</td><td>294.633</td><td>297.9057</td><td>246.4353</td><td>202.4471</td><td>214.702</td></tr><tr><td>4th order polynomial + temperature</td><td>193.0919</td><td>125.6549</td><td>144.8236</td><td>241.0375</td><td>187.4392</td><td>296.628</td><td>299.8974</td><td>247.4832</td><td>198.1059</td><td>215.3708</td></tr><tr><td>5th order polynomial + temperature</td><td>187.4448</td><td>126.6208</td><td>142.6538</td><td>242.7152</td><td>182.5773</td><td>298.1341</td><td>301.8939</td><td>247.9396</td><td>198.5003</td><td>217.049</td></tr><tr><td>6th order polynomial + temperature</td><td>179.8237</td><td>127.0546</td><td>136.7759</td><td>241.5779</td><td>175.0605</td><td><strong>293.6412</strong></td><td>295.469</td><td>249.1004</td><td>198.8159</td><td>215.2838</td></tr><tr><td>Linear model * temperature</td><td>191.4915</td><td>152.4921</td><td>165.7432</td><td>271.1041</td><td>195.8631</td><td>299.7542</td><td>301.6908</td><td>261.5845</td><td>214.1236</td><td>218.4733</td></tr><tr><td>2nd order polynomial * temperature</td><td>191.8085</td><td>140.6064</td><td>144.8647</td><td><strong>238.9733</strong></td><td>186.3509</td><td>300.0211</td><td>301.5847</td><td>249.0721</td><td>216.6558</td><td><strong>208.6308</strong></td></tr><tr><td>3rd order polynomial * temperature</td><td>195.7395</td><td>133.0582</td><td>141.3383</td><td>242.2301</td><td>189.1963</td><td>298.8941</td><td>301.3621</td><td>251.8635</td><td>205.4194</td><td>211.4949</td></tr><tr><td>4th order polynomial * temperature</td><td>191.7946</td><td>112.8464</td><td>140.6139</td><td>241.6344</td><td>173.383</td><td>301.2723</td><td>304.9136</td><td>250.3146</td><td><strong>197.6711</strong></td><td>212.7588</td></tr><tr><td>5th order polynomial * temperature</td><td>188.5724</td><td><strong>112.5853</strong></td><td>138.6162</td><td>243.1846</td><td><strong>170.3905</strong></td><td>300.6816</td><td>304.1199</td><td>252.6209</td><td>198.3532</td><td>214.3907</td></tr><tr><td>6th order polynomial * temperature</td><td>187.5341</td><td>114.5738</td><td>137.9766</td><td>244.8093</td><td>171.4527</td><td>301.4109</td><td>302.7866</td><td>254.381</td><td>200.1839</td><td>214.5539</td></tr></tbody></table>

We used polynomials of various orders (up to 6 since 7 different sampling dates were available) to model gene expression profiles at population scale, accounting for an effect of experimental temperature. We started with the most complex model as follows:

$$ \log(E_{i,j}) = \sum_{i,j} (\beta_{i,j} \times t^i) \qquad (1) $$

where $E_{i,j}$ is the predicted gene expression at the $i^{th}$ sampling date (starting from 0) and at temperature $j$ (considered as a categorical factor). $\beta_{i,j}$ are the parameters of the gene expression model. We then used a stepwise approach to evaluate the fit of less complex models down to the null model (in which $E_{i,j}$ is a scalar) and we computed adjusted R-squared coefficients to compare the amount of explained variation between the selected models for the different genes.

Gene expression dynamics models all exhibited a temporal pattern (Table 1), thereby indicating that none of the considered genes had constant expression through time. Consequently, R² coefficients were all significantly different from 0. Overall, *atg* genes had low R² coefficients (from 18.06% to 46.69%), lysosome-related genes had intermediate to high R² coefficients (30.61% to 86.19%) and metabolism-related genes always had high R² coefficients (73.14% to 87.72%), indicating a stronger inter-individual variation in the former gene groups than in the latter (Fig. 1 and Supplementary online material 1).

Temporal dynamics were always non-linear, with at least a quadratic term (for *ulk1*, *cathl* and *cpt1*) to represent this non-linearity (Table 1, Fig. 1). All other monitored genes displayed a more fluctuating expression pattern with an upward trend for *atg12* and lysosomal genes (*catha*, *cathd* and *cathf*) and a downward trend for *lc3b* and *atg7* (Fig. 1).

Finally, the temperature had no effect on some genes expression (*atg7*, *lc3b*, *cathf* and *cathl*) and interacted with the temporal effect for other genes expression (*atg12*, *ulk1*, *catha*, *gdh* and *cpt1*) (Table 1).

**Inference of transcripts dynamics effects on weight loss.** We hypothesized that gene transcripts quantities can be a limiting factor for cellular functions such as autophagy, especially for long term sustainability⁹,¹⁶,²³⁻²⁵. The usual approach, however, is to estimate simple correlations between punctual levels of gene expression, and the phenotype (here weight loss). We advocate that such approach is unlikely to underpin real relationships between transcripts dynamics and integrated phenotypes, since weight loss is calculated over time (days or weeks of experiment in our case). We therefore aimed at calculating the cumulated production of transcripts over time, based on the gene expression dynamics models previously described and infer the effect of the mean of this quantity on weight loss (shown in Supplementary online material 2).

However, because weight loss is a cumulative data, it is important to notice that any random gene expression profile will be likely to explain a large part of the observed variation in weight loss, especially if i) weight loss is linear over time and ii) the variance of the random gene expression is reduced. The approach used here is to rely on non-linear variation in weight loss, as well as on non-random gene expression profiles. To efficiently assess the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Expression dynamics for autophagy and lysosome related genes. For each gene, blue square symbols indicate data at 9°C, red circle symbols indicate data at 12°C, black triangle symbols indicate data at the beginning of the experiment. Artificial jitter was added on the X coordinates to facilitate data perception. Blue and red lines indicate the median of models predictions at 9°C and 12°C respectively, the pale blue and pale red surfaces represent the 95% confidence interval of the prediction at 9°C and 12°C respectively. When no effect of temperature was detected, the median is indicated by a black line and the surface is displayed in grey. The quality of the model fit for each gene is provided with the adjusted R².</figcaption></figure>

effect of gene expression profiles on weight loss, we therefore generated a null gene expression model depicting the variation of a randomly expressed gene ($E_{null}$), with a variance being the average of observed variances in our dataset. The role of this null gene expression model is to correlate to the linear increase of weight loss, thereby avoiding to mistakenly attribute any linear variation to other genes of interest. The cumulated production of transcripts for the randomly expressed gene at time i and temperature j was expressed as:

$$P(\text{null}, i, j) = \int_{0}^{i} E_{\text{null},i,j} d_i \qquad (2)$$

We then restrained our exploratory analysis to 6 genes of interest: first, we considered 4 mandatory autophagy related genes (ulk1, lc3b, atg7, atg12). If one of these transcripts is lacking in the cell, then the whole construction of the autophagosome is halted. To account for this dependence, we therefore calculated the cumulated production of transcripts for autophagy $P(A)$ at time i and temperature j as the integral of the product of autophagy related transcripts in the cell as follows:

$$P(A, i, j) = \int_{0}^{i} E_{ulk1_{i,j}} \times E_{lc3b_{i,j}} \times E_{atg7_{i,j}} \times E_{atg12_{i,j}} d_i \qquad (3)$$

The calculation of cumulated production of transcripts for *atg* genes ($P(A)$), Supplementary online material 3A) showed a nonlinear accumulation of this production: production was high during the first ten days, then stabilized at a slower rate before raising again during the last ten days. This pattern was more pronounced at high temperature.

Then, we also considered key actors of the catabolism of amino acids and lipids respectively. We calculated the cumulated production of transcripts for *gdh* $P(B)$ and for *cpt1* $P(C)$ at time i and temperature j as follows:

$$P(B, i, j) = \int_{0}^{i} E_{gdh_{i,j}} d_i \qquad (4)$$

$$P(C, i, j) = \int_{0}^{i} E_{cpt1_{i,j}} d_i \qquad (5)$$

For *gdh*, the cumulated production of transcripts ($P(B)$, Supplementary online material 3B) followed a near linear pattern at low temperature, whereas it became nonlinear at high temperature: production was initially higher up to the 20th day, and then paced down to a lower rate. Regarding *cpt1*, the cumulated production of transcripts ($P(C)$, Supplementary online material 3C) displayed a logarithmic pattern, constantly slowing over the duration of the experiment. $P(C)$ was higher at high temperature.
<table><caption>**Table 2.** Relative contributions (reported in %) of *atg* genes, *gdh* and *cpt1* and their interaction in the inference model for weight loss. <sup>a</sup>The contributions have been scaled to 100% after removing the linear trend, that explained 30.64% and 51.71% of the weight loss at 9°C and 12°C respectively.</caption><thead><tr><th rowspan="2">Process</th><th colspan="2">Relative contributions<sup>a</sup></th></tr><tr><th>9°C</th><th>12°C</th></tr></thead><tbody><tr><td>atg genes</td><td>15.31% (SD = 0.137)</td><td>17.84% (SD = 0.149)</td></tr><tr><td>atg genes * gdh</td><td>14.91% (SD = 0.134)</td><td>15.74% (SD = 0.138)</td></tr><tr><td>gdh</td><td>39.33% (SD = 0.19)</td><td>44.6% (SD = 0.205)</td></tr><tr><td>atg genes * cpt1</td><td>11.37% (SD = 0.11)</td><td>6.57% (SD = 0.069)</td></tr><tr><td>cpt1</td><td>19.09% (SD = 0.159)</td><td>15.25% (SD = 0.135)</td></tr></tbody></table>

Then, we assumed that *gdh* and *cpt1* expression could be linked either to autophagy related catabolism, or to autophagy independent catabolism. Indeed, under fasting conditions, protein catabolism (both autophagy-dependent and non-dependent) results in the production of free amino acids, which are used by several enzymes including GDH to fuel the tricarboxylic acid (TCA) cycle for energy production<sup>26</sup>. Similarly, fatty acids produced by lipophagy as well as cytosolic lipases are transported into mitochondria through the carnitine palmitoyltransferase system (including the outer mitochondrial membrane protein CPT1) and enter the fatty acid β-oxidation pathway to produce acetyl-CoA, which fuels the TCA cycle and energy production<sup>7</sup>. The flux of fatty acid β-oxidation pathway is primarily determined by CPT1<sup>27,28</sup>. The general inference of the cumulated production of transcripts on weight loss (W) at time i and temperature j was thus modelled as:

$$W_{i,j} = \gamma_j \times [\alpha_{1,j} \times P(A, i, j) + \alpha_{2,j} \times P(A, i, j) \times P(B, i, j) + \alpha_{3,j} \times P(A, i, j) \\ \times P(C, i, j) + \alpha_{4,j} \times P(B, i, j) + \alpha_{5,j} \times P(C, i, j) + \alpha_{6,j} \times P(null, i, j)] \quad (6)$$

where $\gamma_j$ is a temperature dependent scaling parameter for the general effect of transcripts on weight loss $W_{i,j}$, and $\alpha_{1,j}$, $\alpha_{2,j}$, $\alpha_{3,j}$, $\alpha_{4,j}$, $\alpha_{5,j}$ and $\alpha_{6,j}$ are the relative contributions (i.e., summing to one) of the different transcripts and their interactions sampled in a Dirichlet distribution. $P(C)$ was relatively lower at 12°C than at 9°C (Supplementary online material 3A).

We compared the above presented model (termed as M1) where $P(A)$, $P(B)$, $P(C)$ and $P(null)$ are included, to a model where only the effect of $P(null)$ was inferred on weight loss (M0). When comparing the M0 inference model (where only the random gene expression was used) and the M1 model (where $P(A)$, $P(B)$ and $P(C)$ and their interactions were also used), the DIC value indicated the M1 model as the best compromise between fit and complexity to explain weight loss data ($DIC_{M1} = 921.1$, $DIC_{M0} = 926.4$). In the M1 model, the relative contributions of *atg* genes, *gdh* and *cpt1* on weight loss differed as a function of temperature, although a general pattern emerged. The linear trend explained by randomly expressed gene ($E_{null}$) represented 30.64% of the contribution at 9°C versus 51.71% at 12°C. The general effect of transcript production on weight loss was higher at 12°C ($\gamma = 1.104$, SD = 0.147) than at 9°C ($\gamma = 0.8738$, SD = 0.2038). At 9°C, the biggest contributor to weight loss was $P(B)$ (*gdh*), followed by $P(C)$ (*cpt1*) both independently of $P(A)$. The contribution of *atg* genes $P(A)$ ranked third, closely followed by the interactive effect of $P(A)$ and $P(B)$. The interactive effect of $P(A)$ and $P(C)$ ranked last (Table 2). At 12°C, the relative contribution of $P(B)$, $P(A)$ and their interaction increased, whereas the contribution of $P(C)$, and its interaction with $P(A)$ decreased.

## Discussion

Monitoring the dynamics of biological processes as well as the impact of these dynamics on the establishment of integrated phenotypes is of paramount importance for a better understanding of both the processes themselves and their role at the whole organism-level. In that aim, we conducted a novel approach to (i) build a statistical model for the dynamics of expression of several *atg* and lysosomal genes in glass eel during long-term fasting at two temperatures (9°C and 12°C) and (ii) make use of this dynamic model to test whether or not, and to what extent, the autophagy related transcripts dynamics can contribute to changes in an integrated phenotype – here weight loss.

Our analysis predicts that most of monitored autophagy- and lysosome-related genes display a non-random fluctuating expression pattern during long-term fasting in glass eels. Such a long-term temporal pattern of gene expression dynamics is consistent with the currently accepted idea that gene transcription events control long-term autophagy flux<sup>9,16,23–25</sup>. Indeed, the adaptive autophagic stress response has been referred to as being biphasic. Thus, after a rapid autophagic response, which involves the post-translational modification of cytosolic proteins, a collection of transcription factors upregulates the expression of genes encoding proteins that sustain and regulate autophagic flux. In this regard, the predicted fluctuating expression pattern of autophagy- and lysosomal-related genes in glass eel during long-term fasting could indeed capture some of the successive waves of replenishment of critical proteins of this function.

However, besides the fluctuating pattern, our model also predicts upward or downward expression trends according to the considered genes. It is thus interesting to note that *lc3b*, *atg7* and *cathl* (as well as *ulk1* for the first 40 days of fasting) showed a rather high expression level at the beginning followed by a constant decline, possibly reflecting a fall in autophagic capacity in long-term fasting glass eels. The result obtained for *lc3b* particularly caught our attention. Indeed, previous studies from both yeast and mammalian models supported a positive
correlation between the concentration of LC3b and the mean autophagic vesicle size<sup>29,30</sup>. Notably, mutant yeast strains engineered for reduced *atg8* (LC3b homolog in yeast) expression accumulate autophagosomes smaller than wild-type yeast in response to starvation<sup>30</sup>. Within the mutant strains, the more substantial the *atg8* reduction, the smaller the size of resulting autophagosomes<sup>30</sup>. More recently, Martin *et al.* developed a computational model describing autophagic vesicle dynamics in a mammalian system and also reported a positive correlation between LC3 levels and autophagic vesicle size across single cells, in line with a proposed function for LC3 in phagophore elongation<sup>14</sup>. Collectively, these studies in combination with the temporal downward expression trend of *lc3b* predicted by our model support a decrease in autophagic vesicle size, and therefore in the autophagic capacity, in glass eels during prolonged starvation. These results are consistent with earlier studies showing that protein breakdown through macroautophagy gradually decline in organs such a liver during long-term fasting but remains sustained through Chaperone-Mediated Autophagy (CMA), a selective form of autophagy that does not require the involvement of any of the *atg* genes monitored in the present study<sup>31</sup>. However, it remains to be established whether or not CMA or a CMA-like process exists in fish. In that condition, the predicted upward trend in the expression of some other studied genes (*atg12*, *catha*, *cathd*, *cathf* and also *ulk1* at the end of the experiment) may appear contradictory. However, we and other have previously demonstrated that induction of *atg* and lysosomal transcripts can also occur as a compensatory response to autophagy inhibition<sup>32–35</sup>. The underlying mechanisms are now well described and would at least partially involve the transcription factor TFEB recognized as a major player in the regulation of the expression of *atg* and lysosomal genes<sup>36</sup>. In any case, our analysis highlights the complexity of the regulatory network involved in the autophagy-lysosomal system.

Such a complexity is further revealed when considering the R² coefficient that compares the amount of explained variation between the selected models for the different genes. Indeed, the lower R² coefficients for *atg* genes (from 18.06% to 46.69%) than for lysosomal genes (30.61% to 86.19%) could be explained by the existence of a multiplicity of selective autophagy that converges to a common lysosomal degradative system. Indeed, autophagy was long considered to be a nonselective degradation pathway<sup>37–41</sup>. However, recent research has revealed that autophagy is also critical for the selective degradation of specific cargos, such as aggregated proteins (aggrephagy), endoplasmic reticulum (ER-phagy or Reticulophagy), mitochondria (mitophagy), peroxisomes (pexophagy), lipid droplets (Lipophagy), glycogen (Glycophagy) and even invading pathogens (xenophagy)<sup>42–45</sup>. All these types of selective autophagy require the formation of vesicles and therefore the transcription of the *atg* genes monitored in the present study. However, emerging evidence suggests that the size as well as the number of the autophagosomes not only regulate efficiency, but may also be critical for the cargo selectivity<sup>46</sup>. The relative strong inter-individual variation observed for the monitored *atg* genes may therefore reflect this multiplicity of autophagy-related mechanisms likely activated during fasting.

Finally, contrary to what we might have thought, we do not observe a general trend of the effect of temperature on *atg* and lysosomal gene expression dynamics. The observed effects are relatively slight and concern only a few genes (*ulk1*, *atg12*, *catha* and *cathd*) at certain sampling times. Only one gene (*cathd*) exhibit a slightly higher expression level at 12°C than at 9°C throughout the fasting period. However, these data do not indicate that the dynamics of expression of these genes cannot be affected by temperature, but rather that at the tested temperatures (which correspond to a natural range) few if no effect is observed. Could this apparent low ability of eels to increase the expression of autophagy- and lysosome-related genes under increasing temperature be critical for their migration in the future?

Now, what are the use of these transcripts dynamics when looking at integrated phenotype? We know that autophagy (in its narrow definition of this paper) is only a part of the catabolic chain that provides energy during a fasting period. We therefore wanted to weigh the contribution of the autophagy related transcripts produced over time on an integrated phenotype (here the weight loss) by considering how these genes interact with markers known to drive energy supply through catabolism of proteins and lipids. Looking at the accumulation of *atg* genes, *gdh*, and *cpt1*, it was clear that the dynamics of these 3 groups of transcripts strongly differed during our experiment, which allowed for further inference on weight variation. We found the contribution of autophagy related transcripts dynamics to be significant and quite stable when temperature changed (between 41.58% at 9 °C and 40.15% at 12 °C, summing all terms where autophagy is involved in Table 2). This relatively stable contribution of *atg* genes between the two temperatures may be linked to the wide range of different processes mentioned above (lipophagy, glycophagy, mitophagy, pexophagy etc...) that all depend on the studied genes and that may have been affected differently by temperature. By comparison, the contribution of *cpt1* dramatically reduced at higher temperature, compensated by a relative increase of *gdh* contribution, in line with a preferred use of amino acids relative to fatty acids for energy supply during starvation in European eel<sup>47</sup>. To our knowledge, this is the first demonstration of the effect of *atg* transcripts dynamics on an integrated phenotype. In recent years, a major effort has been made to understand the mechanisms underlying the regulation of *atg* gene expression as well as the role of these nuclear events in the control of cellular homeostasis during long period of stress<sup>9,16,23–25</sup>. However, usual approaches using simple correlations between punctual levels of gene expression and integrated phenotypes failed to capture such functional and dynamic relationships between transcripts and integrated phenotypes.

Overall, the model presented here is consistent with published findings and improved, to some extent, our understanding of weight loss dynamics under long-term fasting and in different temperatures. However, it is currently limited to only minimal molecular details (i.e., 4 transcripts that are mandatory for autophagy). In the future, an extended model, based on the framework presented here and including additional known molecular factors involved in the different types of selective autophagy will certainly provide a better description of temporal individual variation in the integrated phenotype. Additionally, better estimates and results will be obtained using biological models where multiple sampling over time for gene transcription is an option, thereby relaxing our assumption of similarity in the pattern of expression dynamics between individuals, thus simplifying the first stage of our approach.
**Wrapping-up.** Despite these shortcomings, our results hint at the role of *atg* transcripts as potential indicators of ability to cope with long-term fasting in young eels. Linking these results with the sources of variation in gene transcripts dynamics therefore opens exciting research paths to understand (i) finer regulations mechanisms of autophagic response to long-term fasting, (ii) sources of inter-individual variation in ability to cope with environmental variation of resources during migration, a source of growing concern in the various ecosystems occupied by European eels⁴⁸. The adaptive potential of the species is indeed at stake here, since it is also known that determinism of sex is environmentally controlled, especially by growth rate and population densities: on the one hand, the inability of some individuals to sustain long term fasting under changing conditions could affect density distribution along the various ecosystems that stem the European eel life cycle. Yet on the other hand, the control of autophagy on weight loss will interactively affect sex determinism and thereby life cycle length and migration paths. We therefore advocate for further research to focus on long term analysis of gene expression and relationship to phenotypic variation. Based on the present results, we deem it very unlikely that non-random temporal dynamics of gene transcripts, when detected, may contain no information on the dynamics of integrated phenotypes.

**Materials and Methods**

**Origin and handling of fish.** Glass eels were caught in March 2014 (authorization of scientific fishing, 17 September 2013) on the Atlantic coast, downstream from the mouth of the Huchet river (Moliets, South-Western France, 43°51N, 1°23W). Fishing took place by night during flood tide using a hand net. After collection, glass eels were brought to the laboratory and kept in a tank at 11 °C. The following morning, 150 glass eels were anaesthetized by immersion in a solution of benzocaine (0.3 mL L⁻¹), measured (±0.5 mm), weighed (Sartorius CP 153 balance, ±1 mg) and individually tagged using visible implant elastomers as described previously⁴⁹.

**Experimental design.** Ten Aquaria (50 cm long, 25 cm large and 25 cm high, 20 liters) filled up with tap water to 15 cm high (outflow height) were placed in two thermoregulated rooms in order to get five aquaria at 9 °C and the five others at 12 °C. Water was re-circulated through a filter and UV light. Four pieces of artificial vegetation mimicking *Elodae* shelters (10 cm high) were placed on the floor of each aquarium, regularly spaced along a line at 10 cm from the front glass of the aquarium. Photoperiod was 12 L/12D with dim light during the light phase (27 µW/cm² ≈ 4 lux) and 30 minutes of dawn and dusk. Temperature was measured every day and the water level adjusted every week with aerated tap water kept at room temperature.

Fifteen tagged glass eels were transferred in each aquaria (Fourteen for analyses and one additional in case of mortality). There was no initial significant difference between glass eels placed at 9 °C versus 12 °C in weight (Welch 2-samples t-test, t = 1.5198, 131 df, p = 0.131) or length (t-test, t = 0.51773, 131 df, p = 0.6055). Two individuals were randomly sampled in each tank every 10 days during a period of 60 days. Glass eels were immediately anaesthetized, weighted and lengthened as described above, frozen in liquid nitrogen and then kept at -80 °C until analysis. No mortality was observed and the remaining fish were not used for analyses.

Procedures used in this study have been validated by the ethics committee 'Aquitaine poissons-oiseaux', the 15 January 2013. The experiment was carried out in strict accordance with the EU legal frameworks, specifically those relating to the protection of animals used for scientific purposes (i.e., Directive 2010/63/EU), and under the French legislation governing the ethical treatment of animals (Decret no. 2013-118, February 1st, 2013). The INRA/UPPA experimental station is certified for animal services under the permit number B64-495-2 by the French veterinary services, which is the competent authority.

**Behavioral observations.** Fish behavior was observed every morning during the 60 days of experiment. Observations took place at 70 cm away from the front of the aquaria avoiding any visual disturbance. Each aquarium was observed during 5 minutes: Glass eels were identified by their tagging and individual activity (swim or rest) and location (floor, shelter or water column) were noted. At 9 °C, all glass eels remained inactive, hidden in the shelters or on the bottom while at 12 °C, eight glass eels were observed in more of 70% of the observations swimming in the water column or along the wall of the aquarium. As autophagy is related to fasting but also to locomotor activity and stress, these eight glass eels were removed from our analysis.

**mRNA levels analysis: quantitative RT-PCR.** The extraction of total RNA was performed using TRIzol reagent (Invitrogen, 15596018) according to the manufacturer's recommendations. One microgram of the resulting total RNA was reverse transcribed into cDNA, using the SuperScript III RNAsenseH-reverse transcriptase kit (Invitrogen, 18080085) with random primers (Promega, Charbonnières, France) according to the manufacturer's instructions. Primers specific to autophagy- (ulk1, lc3b, atg7, atg12) and lysosome- (cathepsin A, cathepsin D, cathepsin F, cathepsin L) related genes as well as specific to glutamate dehydrogenase (gdh) and carnitine palmitoyltransferase I (cpt1) were designed using Primer3 software (version 4.1.0)⁵⁰,⁵¹. To confirm specificity of the newly developed RT-PCR assay, the amplicon was purified and sequenced (Beckman-Coulter Genomics, Takeley, UK). The primers used for RT-PCR assays are listed in Table 3. Quantitative RT-PCR was carried out on the Roche LightCycler 480 System (Roche Diagnostics, Neuilly sur Seine, France). The assays were performed using a reaction mix of 6 µl per sample, each of which contained 2 µl of diluted cDNA template, 0.24 µl of each primer (10 µM), 3 µl Light Cycler 480 SY Green Master mix (Roche Diagnostics, 4887352001) and 0.52 µl DNAsense/RNAsense-free water (5 Prime GmbH, 2500020). The PCR protocol was initiated at 95 °C for 10 min for initial denaturation of the cDNA and hot-start Taq-polymerase activation, followed by 45 cycles of a 3-step amplification program (15 s at 95 °C; 10 s at 60 to 64 °C and 15 s at 72 °C), according to the primer set used. Melting curves were systematically monitored (temperature gradient at 1.1 °C/10 s from 65 to 94 °C) at the end of the last amplification cycle to confirm the specificity of the amplification reaction. Each PCR assay included replicate samples (duplicate of reverse transcription and PCR amplification, respectively) and negative controls (reverse transcriptase- and
<table><caption>**Table 3.** Sequences of the primer pairs used for real-time quantitative RT-PCR.</caption><thead><tr><th>Gene</th><th>5'/3' Forward primer</th><th>5'/3' Reverse primer</th></tr></thead><tbody><tr><td colspan="3"><strong>Autophagy-related genes</strong></td></tr><tr><td>ulk1</td><td>GACTTTAACGGAGGCGAGTG</td><td>CCAGTCGTGTTTCTCCTTATGTC</td></tr><tr><td>lc3b</td><td>TACAGGACATAGGCCGCTAA</td><td>ACTCGCTGTTCAAATGTCCT</td></tr><tr><td>atg7</td><td>CGGCTGAGATCTGGGACA</td><td>AGCCAGATTGAGCGACTGAT</td></tr><tr><td>atg12</td><td>GCGGTAGGGGACACTCCTAT</td><td>CACTGCCAAAAACATTCAAATAAC</td></tr><tr><td colspan="3"><strong>Lysosome-related genes</strong></td></tr><tr><td>catha</td><td>GGGAACAAGCACCTGCATTA</td><td>CGCCATCATCCTGAATTAGA</td></tr><tr><td>cathd</td><td>GGACATCCCATGCTCTTGTT</td><td>CCGTAATACTGAGCATCAAGGT</td></tr><tr><td>cathf</td><td>GGGATATGGACATCGTAATGG</td><td>GCAGATGGGCGTTGTTTAAT</td></tr><tr><td>cathl</td><td>TCAGTTCTACCAATCTGGAATCTAC</td><td>CTTCCTGTCTTTGGCCATGT</td></tr><tr><td colspan="3"><strong>Metabolism-related genes</strong></td></tr><tr><td>gdh</td><td>TACGCTCAGAACATGGTTGC</td><td>AAGCCGTAGGGGAAAAAGAA</td></tr><tr><td>cpt1</td><td>ATTTGGACCGGTGGCTGAT</td><td>CCTCGAATCACATCTTGTTTCC</td></tr><tr><td colspan="3"><strong>Reference gene</strong></td></tr><tr><td>ef1a1</td><td>CACCGAGACGGAACCTTAAA</td><td>GACCCCTTCCCCATCTCA</td></tr></tbody></table>

cDNA template-free samples, respectively). For the expression analysis of mRNA, relative quantification of target gene expression was performed using the delta CT method described by Pfaffl⁵². The relative gene expression of eukaryotic translation elongation factor 1 alpha 1 (Eef1a1) was used for the normalization of measured mRNA. Its expression showed an acceptable coefficient of variation (CV = 0.296), a positive linear temporal trend was detected but it was too small compared to the temporal variation in expression of other genes to affect our conclusions (Supplementary online material 4). In all cases, PCR efficiency (E) was measured by the slope of a standard curve using serial dilutions of cDNA (1/20, 1/40, 1/80, 1/160, 1/320). In all cases, PCR efficiency values ranged between 1.8 and 2.2.

**Data modelling.** We used R Software (R 3.2.2)⁵³ to fit gene expression polynomials and determinate optimal orders through Akaike Information Criterion (AIC) stepwise selection approach, so to obtain a balance between model fit quality, and complexity of the model (number of parameters). We then modelled gene expression profiles (as deduced from the previous selection) and their effect on weight loss in a single Hierarchical Bayesian model (See Supplementary online material 5), using non informative Gaussian priors for all gene expression models parameters. We used a [0;1000] uniform distribution prior for the $\gamma_i$ parameter, and equiprobable priors (1;1;1;1;1;1) for the Dirichlet distribution. The joint posterior distribution of all hyper parameters of the model was approximated by MCMC sampling as implemented by the OpenBUGS (version 3.2.3) software. A MCMC sample of 10000 draws was used, after checking its convergence by applying the Gelman-Rubin test⁵⁴. The fit of models to observed weight loss data was assessed using the Deviance Information Criterion (DIC)⁵⁵, that balances the quality of the fit (maximum likelihood) by the complexity of the model (number of estimated parameters). Such an approach is best suited when one desires to account for parsimony as much as for predictive power.

**Data availability.** All data generated or analyzed during this study are included in this published article (and its Supplementary Information files).

**References**

1. Dekker, W. Monitoring of glass eel recruitment (RIVO-Netherlands Institute for Fisheries Research, 2002).

2. Tesch, F. The eel: biology and management of anguillid eels (ed. Greenwood, P H.) (John Wiley and Sons, New York, USA, 1977).

3. Daverat, F., Tomas, J., Lahaye, M., Palmer, M. & Elie, P. Tracking continental habitat shifts of eels using otolith Sr/Ca ratios: validation and application to the coastal, estuarine and riverine eels of the Gironde-Garonne-Dordogne watershed. *Marine and Freshwater Research* 56, 619–627, https://doi.org/10.1071/MF04175 (2005).

4. Helfman, G. S., Facey, D. E., Hales, L. S. Jr. & Bozeman, E. L. Jr. Reproductive ecology of the American eel. *American Fisheries Society Symposium*. 1, 42–56 (1987).

5. Oliveira, K. & McCleave, J. Regional variation and the effect of lake: river area on sex distribution of American eels. *Journal of Fish Biology* 58, 943–952 (2001).

6. Bardonnet, A. & Riera, P. Feeding of glass eels (Anguilla anguilla) in the course of their estuarine migration: new insights from stable isotope analysis. *Estuarine Coastal and Shelf Science* 63, 201–209, https://doi.org/10.1016/j.ecss.2004.11.009 (2005).

7. Kaur, J. & Debnath, J. Autophagy at the crossroads of catabolism and anabolism. *Nat Rev Mol Cell Biol* 16, 461–472, https://doi.org/10.1038/nrm4024 (2015).

8. Feng, Y., Yao, Z. & Klionsky, D. J. How to control self-digestion: transcriptional, post-transcriptional, and post-translational regulation of autophagy. *Trends Cell Biol* 25, 354–363, https://doi.org/10.1016/j.tcb.2015.02.002 (2015).

9. Fullgrabe, J., Ghislat, G., Cho, D. H. & Rubinsztein, D. C. Transcriptional regulation of mammalian autophagy at a glance. *J Cell Sci* 129, 3059–3066, https://doi.org/10.1242/jcs.188920 (2016).

10. Cuervo, A. M. Cell biology. Autophagy's top chef. *Science* 332, 1392–1393, https://doi.org/10.1126/science.1208607 (2011).

11. Borlin, C. S., Lang, V., Hamacher-Brady, A. & Brady, N. R. Agent-based modeling of autophagy reveals emergent regulatory behavior of spatio-temporal autophagy dynamics. *Cell Commun Signal* 12, 56, https://doi.org/10.1186/s12964-014-0056-8 (2014).

12. Holczer, M., Marton, M., Kurucz, A., Banhegyi, G. & Kapuy, O. A Comprehensive Systems Biological Study of Autophagy-Apoptosis Crosstalk during Endoplasmic Reticulum Stress. *Biomed Res Int* 2015, 319589, https://doi.org/10.1155/2015/319589 (2015).
13. Kapuy, O., Vinod, P. K., Mandl, J. & Banhegyi, G. A cellular stress-directed bistable switch controls the crosstalk between autophagy and apoptosis. *Mol Biosyst* **9**, 296–306, https://doi.org/10.1039/c2mb25261a (2013).

14. Martin, K. R. et al. Computational model for autophagic vesicle dynamics in single cells. *Autophagy* **9**, 74–92, https://doi.org/10.4161/auto.22532 (2013).

15. Tavassoly, I. et al. Dynamic Modeling of the Interaction Between Autophagy and Apoptosis in Mammalian Cells. *CPT Pharmacometrics Syst Pharmacol* **4**, 263–272, https://doi.org/10.1002/psp4.29 (2015).

16. Lapierre, L. R., Kumsta, C., Sandri, M., Ballabio, A. & Hansen, M. Transcriptional and epigenetic regulation of autophagy in aging. *Autophagy* **11**, 867–880, https://doi.org/10.1080/15548627.2015.1034410 (2015).

17. Xie, Z. & Klionsky, D. J. Autophagosome formation: core machinery and adaptations. *Nat Cell Biol* **9**, 1102–1109, https://doi.org/10.1038/ncb1007-1102 (2007).

18. Punnonen, E. L., Autio, S., Marjomaki, V. S. & Reunanen, H. Autophagy, cathepsin L transport, and acidification in cultured rat fibroblasts. *J Histochem Cytochem* **40**, 1579–1587, https://doi.org/10.1177/40.10.1326577 (1992).

19. Turk, V., Turk, B., Guncar, G., Turk, D. & Kos, J. Lysosomal cathepsins: structure, role in antigen processing and presentation, and cancer. *Adv Enzyme Regul* **42**, 285–303, https://doi.org/10.1016/S0065-2571(01)00034-6 (2002).

20. Uchiyama, Y. Autophagic cell death and its execution by lysosomal cathepsins. *Arch Histol Cytol* **64**, 233–246 (2001).

21. Mommsen, T. P. Salmon spawning migration and muscle protein metabolism: the August Krogh principle at work. *Comp Biochem Physiol B Biochem Mol Biol* **139**, 383–400, https://doi.org/10.1016/j.cbpc.2004.09.018 (2004).

22. IPCC, 2014: Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Core Writing Team, R. K. Pachauri & L. A. Meyer (eds)]. IPCC, Geneva, Switzerland, 151 pp.

23. Botti-Millet, J., Nascimbeni, A. C., Dupont, N., Morel, E. & Codogno, P. Fine-tuning autophagy: from transcriptional to posttranslational regulation. *Am J Physiol Cell Physiol* **311**, C351–362, https://doi.org/10.1152/ajpcell.00129.2016 (2016).

24. Fullgrabe, J., Klionsky, D. J. & Joseph, B. The return of the nucleus: transcriptional and epigenetic control of autophagy. *Nat Rev Mol Cell Biol* **15**, 65–74, https://doi.org/10.1038/nrm3716 (2014).

25. Pietrocola, F. et al. Regulation of autophagy by stress-responsive transcription factors. *Semin Cancer Biol* **23**, 310–322, https://doi.org/10.1016/j.semcancer.2013.05.008 (2013).

26. Plaitakis, A., Kalef-Ezra, E., Kotzamani, D., Zaganas, I. & Spanaki, C. The Glutamate Dehydrogenase Pathway and Its Roles in Cell and Tissue Biology in Health and Disease. *Biology (Basel)* **6**, https://doi.org/10.3390/biology6010011 (2017).

27. Eaton, S. Control of mitochondrial beta-oxidation flux. *Prog Lipid Res* **41**, 197–239 (2002).

28. McGarry, J. D. & Brown, N. F. The mitochondrial carnitine palmitoyltransferase system. From concept to molecular analysis. *Eur J Biochem* **244**, 1–14 (1997).

29. Nakagawa, I. et al. Autophagy defends cells against invading group A Streptococcus. *Science* **306**, 1037–1040, https://doi.org/10.1126/science.1103966 (2004).

30. Xie, Z., Nair, U. & Klionsky, D. J. Atg8 controls phagophore expansion during autophagosome formation. *Mol Biol Cell* **19**, 3290–3298, https://doi.org/10.1091/mbc.E07-12-1292 (2008).

31. Tasset, I. & Cuervo, A. M. Role of chaperone-mediated autophagy in metabolism. *FEBS J* **283**, 2403–2413, https://doi.org/10.1111/febs.13677 (2016).

32. Emanuel, R. et al. Induction of lysosomal biogenesis in atherosclerotic macrophages can rescue lipid-induced lysosomal dysfunction and downstream sequelae. *Arterioscler Thromb Vasc Biol* **34**, 1942–1952, https://doi.org/10.1161/ATVBAHA.114.303342 (2014).

33. Pi, H. et al. Enhancing lysosomal biogenesis and autophagic flux by activating the transcription factor EB protects against cadmium-induced neurotoxicity. *Sci Rep* **7**, 43466, https://doi.org/10.1038/srep43466 (2017).

34. Seiliez, I. et al. Looking at the metabolic consequences of the colchicine-based *in vivo* autophagic flux assay. *Autophagy* **12**, 343–356, https://doi.org/10.1080/15548627.2015.1117732 (2016).

35. Settembre, C. et al. A lysosome-to-nucleus signalling mechanism senses and regulates the lysosome via mTOR and TFEB. *EMBO J* **31**, 1095–1108, https://doi.org/10.1038/emboj.2012.32 (2012).

36. Sardiello, M. et al. A gene network regulating lysosomal biogenesis and function. *Science* **325**, 473–477, https://doi.org/10.1126/science.1174447 (2009).

37. Ashford, T. P. & Porter, K. R. Cytoplasmic components in hepatic cell lysosomes. *J Cell Biol* **12**, 198–202 (1962).

38. De Duve, C. & Wattiaux, R. Functions of lysosomes. *Annu Rev Physiol* **28**, 435–492, https://doi.org/10.1146/annurev.ph.28.030166.002251 (1966).

39. Mortimore, G. E. & Schworer, C. M. Induction of autophagy by amino-acid deprivation in perfused rat liver. *Nature* **270**, 174–176 (1977).

40. Novikoff, A. B. & Essner, E. Cytolysomes and mitochondrial degeneration. *J Cell Biol* **15**, 140–146 (1962).

41. Takeshige, K., Baba, M., Tsuboi, S., Noda, T. & Ohsumi, Y. Autophagy in yeast demonstrated with proteinase-deficient mutants and conditions for its induction. *J Cell Biol* **119**, 301–311 (1992).

42. Anding, A. L. & Baehrecke, E. H. Cleaning House: Selective Autophagy of Organelles. *Dev Cell* **41**, 10–22, https://doi.org/10.1016/j.devcel.2017.02.016 (2017).

43. Delbridge, L. M., Mellor, K. M., Taylor, D. J. & Gottlieb, R. A. Myocardial autophagic energy stress responses—macroautophagy, mitophagy, and glycophagy. *Am J Physiol Heart Circ Physiol* **308**, H1194–1204, https://doi.org/10.1152/ajpheart.00002.2015 (2015).

44. Liu, K. & Czaja, M. J. Regulation of lipid stores and metabolism by lipophagy. *Cell Death Differ* **20**, 3–11, https://doi.org/10.1038/cdd.2012.63 (2013).

45. Nakatogawa, H. & Mochida, K. Reticulophagy and nucleophagy: New findings and unsolved issues. *Autophagy* **11**, 2377–2378, https://doi.org/10.1080/15548627.2015.1106665 (2015).

46. Jin, M. & Klionsky, D. J. Regulation of autophagy: modulation of the size and number of autophagosomes. *FEBS Lett* **588**, 2457–2463, https://doi.org/10.1016/j.fecslet.2014.06.015 (2014).

47. Lovern, J. A. Captive eels: some observations on their behaviour. *Salm. Trout Mag.* **94**, 56–57 (1939).

48. Miller, M. J. & Tsukamoto, K. The ecology of oceanic dispersal and survival of anguillid leptocephali. *Canadian Journal of Fisheries and Aquatic Sciences* **74**, 958–971, https://doi.org/10.1139/cjfas-2016-0281 (2017).

49. Imbert, H., Beaulaton, L., Rigaud, C. & Elie, P. Evaluation of visible implant elastomer as a method for tagging small European eels. *Journal of Fish Biology* **71**, 1546–1554, https://doi.org/10.1111/j.1095-8649.2007.01617.x (2007).

50. Koressaar, T. & Remm, M. Enhancements and modifications of primer design program Primer3. *Bioinformatics* **23**, 1289–1291, https://doi.org/10.1093/bioinformatics/btm091 (2007).

51. Untergasser, A. et al. Primer3 - new capabilities and interfaces. *Nucleic Acids Res* **40**, e115, https://doi.org/10.1093/nar/gks596 (2012).

52. Pfaffl, M. W., Horgan, G. W. & Dempfle, L. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. *Nucleic Acids Res* **30**, e36 (2002).

53. R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2016).

54. Gelman, A. & Rubin, D. B. Inference from Iterative Simulation Using Multiple Sequences. *Statist. Sci.* **7**, 457–472, https://doi.org/10.1214/ss/1177011136 (1992).

55. Spiegelhalter, D. J., Best, N. G., Carlin, B. R. & van der Linde, A. Bayesian measures of model complexity and fit. *Journal of the Royal Statistical Society Series B-Statistical Methodology* **64**, 583–616, https://doi.org/10.1111/1467-9868.00353 (2002).
**Acknowledgements**

This study was supported by the University of Pau et Pays de l'Adour (BQR MIRA "Etude de la mise en place des mécanismes de protéolyse chez la civelle d'anguille en période de jeûne. Relation avec leur comportement de migration estuarienne"). We wish to thank Pascale Coste, Emmanuel Huchet and Jacques Rives (from INRA-UPPA UMR 1224 Ecobiop, St Pée sur Nivelle) and Elisabeth Plagnes-Juan (from INRA-UPPA UMR 1419 NuMeA, Saint Pée sur Nivelle) for their technical assistance.

**Author Contributions**

V.B. and I.S. designed and managed the study. L.O. performed behavioral observations as well as gene expression analyses. J.L. performed data modeling. I.S. and J.L. wrote the manuscript. S.P. and V.B. contributed to the manuscript correction.

**Additional Information**

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-017-18164-6.

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017