BACKGROUND
the ability to massively downregulate metabolic and physiological functions during extended periods of fasting has evolved in multiple species of snakes.
this downregulation of physiological form includes the atrophy of organs such as the heart, kidney, liver, and small intestine.
upon feeding, the size and function of these organs, along with oxidative metabolism, is massively upregulated to accommodate digestion .
of the snake species that experience these large fluctuations in physiology, the burmese python  is the most well-studied  <cit> .
within 48 h of feeding, burmese pythons can undergo up to a 44-fold increase in metabolic rate and >100-fold increases in plasma triglyceride content  <cit> .
major organs also experience dramatic shifts in physiological form, including 40–100% increases in the mass of the heart, liver, pancreas, kidneys, and small intestine .
this extreme organ regenerative growth and atrophy is unparalleled across vertebrates, and studies indicate that this organ growth is driven by multiple cellular processes, including cellular hypertrophy in the heart and mixtures of hypertrophy and hyperplasia in the kidney, liver, and small intestine  <cit> .
organ growth peaks around 1– <dig> days post-feeding , and by 10–14dpf, organ form and function, as well as gene expression patterns, have completely reversed back to fasted levels .
previous studies have examined aspects of this post-feeding response using morphological and physiological assays , analyses of gene expression  <cit> , and combinations of the two  <cit> .
together, these studies have demonstrated that transcriptional responses following feeding are extremely rapid and massive, both in the magnitude of expression changes and in the number of genes with significant differential expression.
genes important in a number of developmental, metabolic, proliferative, apoptotic, and growth processes have been shown to be involved in these major shifts in organ form and function  <cit> .
previous studies have shown that mammalian cells respond to the growth signals in post-fed python serum, which likely indicates a conserved response to core signaling molecules  <cit> .
we therefore hypothesize that a relatively small number of core molecular regulatory molecules and signaling pathways may underlie these responses.
however, the identification of a core set of upstream regulatory molecules and mechanisms has been hindered by the large number of genes that are significantly differentially expressed during this response, making manual interpretation of this gene expression data difficult.
additionally, the lack of comparable replicated sampling across multiple organs has further prevented meaningful across-organ comparisons of changes in gene expression in previous studies  <cit> .
accordingly, major gaps remain in our understanding of the specific mechanisms and growth pathways that are responsible for driving these extreme shifts in burmese python organ size and function, as well as how these mechanisms may vary across different organ systems.
our previous study of the burmese python feeding response addressed some gaps through the use of increased replicates and more frequent time point sampling for one organ, the small intestine  <cit> .
we identified over  <dig>  genes that were significantly differentially expressed during post-feeding regeneration in the small intestine with many of these genes being functionally linked to cellular processes such as wnt signaling, cell cycling, and apoptosis.
this study also linked changes in gene expression with functional and phenotypic shifts by comparing rnaseq data with physiological and histological data.
this detailed analysis was only conducted on the small intestine, however, and failed to address any upper-level signaling mechanisms and pathways.
here, we leverage fully replicated organ-specific time courses detailing gene-level responses to infer canonical pathways and regulatory molecules driving post-feeding organ growth in the burmese python.
we examined gene expression across four major organ systems – the heart, liver, kidney, and small intestine.
we combined increased replicated sampling with statistical inferences of pathway activation and regulatory molecule prediction to identify the mechanistic drivers of cross-tissue, post-feeding organ regeneration.
despite highly organ-specific gene expression responses associated with organ regenerative growth, we found evidence for high degrees of overlap in predicted pathways and regulatory molecules underlying these growth processes between organs.
pathways predicted to be involved in regulating this physiological response include lxr/rxr activation, pi3k/akt, and mtor signaling.
interestingly, we also found strong and consistent evidence for the involvement of nrf2-mediated oxidative stress response and other stress-response pathways in this extreme example of rapid organ growth.
our results suggest that post-feeding, regenerative organ growth in the burmese python may stem from small numbers of key effector molecules mediating a core set of growth and stress-response pathways, which in turn activate diverse, tissue-specific signaling cascades.
methods
feeding experiments
burmese pythons were obtained from commercial breeders.
all animal care and tissue sampling was conducted using protocols approved by the university of alabama institutional animal care and use committee .
burmese pythons were sampled at three physiological states: fasted , 1 day post-feeding  and 4dpf, with the meal consumed equaling at least 25% of their body mass.
previous studies have shown that organ masses and functional phenotypes climax between  <dig> and  <dig> dpf  <cit>  and that phenotypes begin to decline by 4dpf  <cit> .
we therefore chose sampling time points here to capture gene expression patters during the period before phenotypes climax  and early in their regression .
snakes were humanely euthanized by severing the spinal cord immediately behind the head, and organs were immediately extracted, snap frozen in liquid nitrogen, and stored at -80 °c.
between three and six biological replicates  were sampled for each time point.
see additional file  <dig>  supplementary methods for additional details.
transcriptome library generation
total rna was extracted from ~50 mg of snap-frozen tissue using trizol reagent , followed by mechanical cell disruption using a tissuelyzer for 10 min at  <dig> strokes/min, and precipitation of rna using isopropanol.
individual illumina mrnaseq libraries were constructed using either the illumina truseq rnaseq kit or the neb next rnaseq kit, both of which included poly-a selection, rna fragmentation, cdna synthesis, and indexed illumina adapter ligation.
completed rnaseq libraries were quantified on a bioanalyzer , pooled in equal molar ratios in various multiplex arrangements, and sequenced on either an illumina gaiix or illumina hiseq <dig> .
quantifying and visualizing gene expression
raw demultiplexed illumina rnaseq reads were quality filtered and trimmed with trimmomatic v.  <dig>   <cit> .
in instances where the same library was sequenced in multiple different runs, reads were combined and mapped for each individual and time point.
mapping of reads to the reference transcriptome of the burmese python  <cit>  was conducted using bwa v.  <dig> .
<dig>  <cit>  with the following parameters: mismatch penalty =  <dig>  gap open penalty =  <dig>  and alignment score minimum =  <dig>  expression was determined using samtools v.  <dig> .
<dig>  <cit>  by counting the number of unique gene reads that mapped to an annotated transcript, while excluding reads that mapped to multiple positions.
new rnaseq data for various time points and replicates was analyzed together with previously published data from other individuals and replicates  <cit> .
newly-generated sequencing data were archived on the ncbi short read archive .
raw expression counts were normalized using tmm normalization in edger  <cit>  and all statistical analyses of gene expression were conducted using normalized data.
we identified genes that were significantly differentially expressed between time points using two approaches.
first, we estimated significant changes in gene expression between pairs of time points using pairwise exact tests for the binomial distribution calculated in edger, integrating both common and tagwise dispersion  <cit> .
second, to accommodate the time-series nature of the experimental design, we also conducted step-wise regression analysis of gene expression in masigpro  <cit> .
regression analysis enabled the detection of genes with significant patterns of differential expression across all three time points.
gene expression heatmaps were generated in r and clustered with the package vegan  <cit> , with gene clustering calculated using average linkage hierarchical clustering based on a bray-curtis dissimilarity matrix.
we used the program stem  <cit>  to identify and visualize significant expression profiles for all genes in our rnaseq data.
assigning homology for functional analyses
to facilitate the use of various pathway activation and regulatory molecule predictions, we annotated the full burmese python transcript set  <cit>  with orthologous human gene ensembl  <cit>  identifiers.
reciprocal tblastx was first conducted between anolis carolinensis and burmese python, and anolis gene ids identified as orthologous to python genes were converted to human ensembl identifiers using homology tables from ensembl’s biomart  <cit> .
the same process of reciprocal best blast using tblastx was performed between burmese python and gallus gallus, followed by conversion of chicken ensembl identifiers to human ensembl identifiers using homology tables from ensembl’s biomart  <cit> .
we also performed reciprocal best blast of the python with homo sapiens.
finally, we used one-way tblastx with anolis, chicken, and human to annotate python genes that were not assigned an ortholog from reciprocal best blast.
using this annotation approach, we were able to assign human ensembl ids to  <dig>  of  <dig>  total python reference transcripts.
pathway and upstream regulatory molecule analyses
to infer the involvement of upstream regulatory molecules and pathways, we performed core analysis in ingenuity pathway analysis , using default parameters.
ipa uses gene identifiers and the fold-change value for each differentially expressed gene to identify enrichment patterns for canonical pathway analysis  and upstream regulatory molecule analysis , and to infer the activation direction  between particular time points.
these two analyses both use observed gene expression data to infer unobserved features , but differ fundamentally in how they use expression data to make inferences.
cpa predicts the involvement and activation/inhibition of canonical pathways based on observed evidence from gene expression data, specifically for genes that participate as higher-level regulatory molecules within a given pathway; analysis of observed gene expression data incorporates information from the ingenuity knowledge base  to provide both a statistical value of enrichment and a prediction of the biological involvement for the pathway as a whole .
in contrast, urma uses observed changes in gene expression specifically for genes at lower levels within pathways  to predict activation or inhibition of regulatory molecules upstream of these genes  <cit> .
due to differences in these approaches, together these two methods provide a well-rounded set of comparable inferences for dissecting molecular mechanisms .fig.
<dig> conceptual overview of differences between canonical pathway analysis  and upstream regulatory molecule analysis .
pairwise analyses on experimental gene expression data  identify significantly upregulated and downregulated genes .
significantly differentially expressed genes are then analyzed in two distinct ipa analyses   canonical pathway analysis predicts pathway activation based on overlap of gene expression data with molecules within the pathway.
d upstream regulatory molecule analysis predicts activation of specific regulatory molecules based on downstream molecules in our gene expression dataset
for ipa analyses, we used only genes identified as significant in pairwise differential expression analyses between time intervals , and we input fold changes per gene averaged across biological replicates, along with our estimate of the orthologous human ensembl id for each gene.
pathways important to cross-tissue physiological responses were isolated using the ipa cpa , with a right-tailed fisher’s exact test p-value of less than  <dig> .
we examined only those pathways that were significant, based on a predicted activation z-score, in at least one of the four organs for at least one of the post-feeding time points.
for ipa analyses, the z-score is used to determine the statistical significance of the number of activated and inhibited predictions, and the sign of the value indicates the overall activation state .
we used a p-value cutoff of  <dig>  for the cpa in ipa to reduce potentially spurious inferences.
upstream regulators and hypotheses for global signaling molecules were identified using urma in ipa, with a fisher’s exact test overlap p-value threshold of  <dig> .
pathway network figures were modified manually from predicted network figures generated in ipa.
for analysis of specific pathways , we also determined the number of genes involved in each pathway that were assigned python orthologs by our orthology analyses, and how many of these genes were expressed at some level in our dataset .
RESULTS
trends in gene expression across organs
we used our expression data from all python samples  to examine the degree to which different organ systems ‘turn on’ upon feeding and then experience ‘regression’ towards pre-feeding patterns of expression at 4dpf.
we found that for each organ, the majority of differentially expressed genes showed immediate up- or downregulation from fasting to 1dpf.
interestingly, each of the four organs examined appeared to experience regression towards fasting levels of expression by 4dpf to widely different extents, indicating that each organ may have its own unique temporal program of growth followed by atrophy.
across organs, the heart appeared to shift towards regression the fastest.
other organs experienced reversals of fasted to 1dpf expression shifts to varying degrees by 4dpf, ranging from the moderately paced small intestine and kidney, to the slow-paced liver .
stem analysis further supported these temporal patterns of up-regulation and regression across organs .table  <dig> numbers of differentially expressed genes between pre- and post-feeding time points for the four organs studied
for each comparison, the numbers of up and downregulated genes were inferred using pairwise analysis with a benjamini-hochberg corrected p-value < <dig>
regression analysis across time points, which tends to be conservative, identified hundreds of genes that were significantly differentially expressed across all three time points with  <dig> genes in the heart,  <dig> genes in the kidney,  <dig> genes in the liver, and  <dig>  genes in the small intestine.
of the  <dig>  total genes differentially expressed across all four organs, 21% are unique to the heart, 16% are unique to the kidney, 15% are unique to the liver, and 32% are unique to the small intestine .
only a single gene was identified as significant in all four organs across all time points: coagulation factor x .fig.
<dig> summary of significantly differentially expressed genes for all four organs identified via regression analysis.
a venn diagram depicting the numbers of genes significantly differentially expressed across time points.
darker colors indicate a large number of genes and lighter colors indicate a smaller number of genes.
b heatmaps depicting all significantly differentially expressed genes across all time points in each organ.
<dig> genes were significantly differentially expressed in the heart.
there were  <dig> genes significantly differentially expressed in the kidney.
<dig> genes were significantly differentially expressed in the liver and  <dig>  genes showed significant differential expression in the small intestine
to further dissect patterns of expression change following feeding, we conducted pairwise analyses of gene expression between time points for each organ.
in the heart, pairwise analyses identified  <dig> significantly differentially expressed genes between fasted and 1dpf , and  <dig> genes were significantly differentially expressed between 1dpf and 4dpf .
in the kidney,  <dig> genes were significantly differentially expressed between the fasted state and 1dpf , while only eight genes were significantly differentially expressed from 1dpf to 4dpf .
in contrast to the heart, we found many genes  significantly differentially expressed between fasted and 4dpf in the kidney.
in the liver,  <dig> genes were differentially expressed within 1dpf , while only  <dig> genes were significantly differentially expressed from 1dpf to 4dpf .
with  <dig> genes significantly differentially expressed between fasted and 4dpf, among all four organs, the liver was the least ‘reset’ to the fasting condition by 4dpf.
finally, the small intestine showed higher levels of differential expression than the other three organs.
within 1dpf,  <dig>  genes were significantly differentially expressed .
from 1dpf to 4dpf,  <dig> genes were upregulated and  <dig> genes were downregulated, and  <dig> genes were differentially expressed between fasted and 4dpf .
genes and pathways implicated in differential gene expression in individual tissues
to move beyond gene-specific responses and towards deciphering the mechanisms that may underlie growth responses across different organs, we identified pathways that were significantly activated/repressed between fasting and 1dpf .
we found consistent evidence that the nrf <dig> stress-response pathway is activated in all tissues, except in the heart, where there was insufficient data to determine the direction of activation.
we also found relatively consistent evidence for activation of the related growth pathways mtor and pi3k/akt across organs, although this inference was most significant in the heart and small intestine.
we also inferred the involvement of the related pathways: lxr/rxr, lps/il-1-mediated inhibition of rxr function, ppar/rxr, and ppar signaling in multiple organs; the direction of stimulation of these pathways was both variable across organs and inconclusive in some organs.
substantial involvement of cytoskeletal pathways, including actin cytoskeleton signaling and actin nucleation by arp-wasp complex, was also inferred across organs and positive in the kidney and small intestine, yet negative or inconclusive for the heart and liver, respectively.fig.
<dig> canonical pathways predicted to be activated or inhibited from gene expression data.
each pathway shown is significantly enriched for our genes with a fisher’s exact test p-value less than  <dig>  .
pathways were shown only if they met our criteria for significance and had a predicted activation state in at least one organ.
z-scores of  <dig>  indicate pathway predictions that lack a bias in the direction of gene regulation observed in our dataset.
ppar signaling  was also included
in addition to pathway activation/repression patterns shared across organs, a number of pathways showed substantial organ-specific directionality of response.
examples of this pattern include the growth-related ampk signaling pathway , erk <dig> signaling , and integrin signaling .
lastly, a number of pathways appeared to be organ-specific, including p <dig> mapk and erk <dig> signaling in the heart and 14-3-3-mediated signaling in the small intestine .
upstream regulatory molecule analysis of  <dig> dpf responses
our inferences of upstream regulatory molecules  between the fasted and  <dig> dpf time points supported many of the same molecular mechanisms underlying organ growth identified via cpa, such as stress response, growth, and lipid signaling pathways.
we explored urm predictions for all classes of urms except biological drugs, chemicals, and micrornas.
we found that many predicted urms were shared among organs, with  <dig> shared among all four organs.
predicted urms also showed substantial organ-specific patterns, with a large number of urms uniquely predicted for each organ.
the heart showed the largest number of unique urm predictions , while only  <dig>   <dig>  and  <dig> unique urms were predicted in the kidney, liver, and small intestine, respectively .fig.
<dig> predicted upstream regulators from ipa analysis of gene expression changes from fasted to 1dpf.
a venn diagram of all upstream regulatory molecules analyzed.
b heatmap of predicted activation z-scores for selected classes of upstream regulatory molecules.
green indicates predicted activation, blue indicates predicted inhibition, white indicates the regulator is not predicted to function in that organ, and grey indicates that the upstream regulator is predicted to have significant involvement but the activation state cannot be determined based on the gene expression data.
regulators shown in this heatmap were filtered by three conditions: 1) were present in at least three of the four organs, 2) are significantly predicted , and 3) have activation z-scores greater than | <dig> | in at least one organ.
biological drug, chemical, and microrna categories were excluded from urm analyses
to identify regulators with broadly relevant patterns across multiple organs, we focused on urms predicted significantly in at least three organs and with moderate to high activation z-score  in at least one organ.
a subset of the urms meeting these criteria is shown in fig.
4b, and the full set is shown in additional file 1: figure s <dig>  many of these urm predictions coincided directly with predicted canonical pathways.
nfe2l <dig> and atf <dig>  key regulators within the nrf2-mediated oxidative stress response pathway, were predicted to be strongly activated in the small intestine, liver, and kidney, consistent with the canonical pathway analysis predictions of activation of the overall nrf <dig> pathway in these three organs.
we also predicted involvement of nfkb and nfkbia, two key regulators within the nfkb signaling response pathway – this inflammatory response pathway is thought to be inhibited by activation of the nrf2-mediated oxidative stress response pathway  <cit> .
nfkb was predicted to be inhibited in the liver and heart, weakly activated in the kidney, and absent in the small intestine, while nfkbia was predicted to be inhibited in the liver, weakly activated in the heart and kidney, and again absent in the small intestine.
activation of the growth pathways mtor and pi3k/akt were additionally supported by activation of predicted regulators such as mtorc <dig> and raf <dig>  respectively, and the inhibition of pten.
lipid signaling pathways such as lxr/rxr signaling, lps/il-1-mediated inhibition of rxr function, ppar/rxr, and ppar signaling were supported by several predicted urms such as pxr ligand, nr1h <dig>  nr1i <dig>  nr1i <dig>  srebf <dig>  srebf <dig>  ppara, pparg, rxra, ppargc1a, and ppargc1b .
these urms were consistently predicted as activated in the small intestine, liver, and kidney and either absent or predicted as inhibited in the heart.
it is notable that while the inferences of activation directions of lipid signaling pathways across organs were largely ambiguous and sometimes inconsistent in our cpa , the associated urms display a consistent trend of predicted activation in the small intestine, liver, and kidney, and either predicted inhibition or absence in the heart .
additionally, several urms, particularly for the mtor pathway, were predicted as inconsistent or even contradictory to the results of cpa or our experimental data.
for example, while mtorc <dig> is predicted as significantly activated in the small intestine by urma, this molecule is downregulated in our experimental data .
additionally, the mtor protein that is involved in forming both of the main complexes of the mtor signaling pathway  is predicted to be strongly inhibited in the small intestine and weakly inhibited in the kidney and heart.
both of these urm predictions appear to contradict the positive activation of the mtor signaling pathway inferred for the small intestine and heart as inferred from the cpa.
in addition to urms involved in key predicted canonical pathways, upstream regulatory analysis predicted several other notable urms with strong activation or informative trends across organs.
insulin and insr were both predicted as strongly activated regulators in the kidney, liver, and small intestine, suggesting a possible role of insulin receptor signaling in facilitating this regenerative response, which is also consistent with activation of the mtor pathway.
myc, a regulator within the erk <dig> and p <dig> mapk signaling pathways, was predicted as activated in all four organs, although strongest in the liver.
several regulators within the mapk signaling pathway were also predicted in urma, with atf <dig> and atf <dig> predicted as activated in the kidney, liver, and small intestine, and erk predicted as activated in the kidney and inhibited in the heart and liver.
these urms suggest the involvement of the erk and mapk signaling pathways in this response, even though cpa predictions for these two pathways were not substantially strong .
detailed dissection of nrf <dig> and mtor pathway responses to feeding
we were particularly interested in our findings that the nrf <dig> stress response and the mtor growth pathways appear to be involved in post-feeding growth in multiple organs.
to investigate these inferences further, we fully dissected evidence from our gene expression data for activation of these pathways by visualizing observed and inferred evidence for activation of these pathways in the context of ipa-generated pathway maps .
specifically, we generated pathway predictions that integrate both observed shifts in gene expression from our data , and estimates of activation/inhibition of molecules downstream of these observed genes that are inferred based on canonical signaling patterns in these pathways.
relevant to our power to detect pathway-wide signals of activity, we were able to associate over 70% of human genes within the mtor and nrf <dig> pathways with python orthologs that were expressed at some level in our dataset ; thus, we expect that our power and degree of resolution of pathway activation for these particular pathways is quite good.fig.
<dig> combined gene expression and predicted activation information for the mtor pathway in the heart and small intestine.
a gene expression and predicted activity for the mtor pathway in the heart.
b gene expression and predicted activity for the mtor pathway in the small intestine.
differentially expressed genes identified in our rnaseq data set are highlighted in red  and blue  while predicted activation states are highlighted in orange  and green .
c cpa and urma results for pathways and upstream regulatory molecules involved in mtor signaling and other relevant growth pathways
fig.
<dig> ipa generated pathway prediction for the nrf2-mediated oxidative stress response in the small intestine.
predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set.
a gene expression and predicted activity for the nrf﻿ <dig> pathway in the small intestine.
differentially expressed genes identified in our rnaseq data set are highlighted in red  and blue  while predicted activation states are highlighted in orange  and green .
b cpa and urma results for pathways and upstream regulatory molecules involved in nrf <dig> signaling and other related pathways
pathways maps for mtor responsiveness between fasted and 1dpf show both common and divergent patterns of pathway activation among organs .
the heart  and kidney  both show similar patterns of mtor activation, including the activation of both the mtorc <dig> and mtorc <dig> complexes.
major differences in mtor activation between these two organs includes strong evidence for downregulation of ampk and the eif <dig> complex in the heart, yet, no direct and/or clear evidence for up- or downregulation of these complexes in the kidney.
in the small intestine, the mtor pathway was inferred to be strongly downregulated, as is akt; ampk and the eif <dig> complex showed mixed signs of activation  .
it is also notable that different organs showed different levels of internal consistency in the integration of results with the known functionality within the mtor pathway.
for example, the heart and kidney have either zero or one pathway connection in which gene expression results contradict the direction of activation of the pathway  – for the kidney this disagreement occurs in the relationship between rsk and inhibition of tsc <dig> .
in the small intestine, eight such disagreements occur , and most of these occur at the steps immediately above or below activation of mtorc <dig> and mtorc <dig> complexes.
the liver was the only organ that contained no signal for the activation or repression of mtor pathway .
it should be noted that inferences for mtor activation from cpa are at times contradictory to those identified via urma .
while predictions based on the pathway maps indicate downregulation of mtor in the small intestine, the z-score suggests slight upregulation of this pathway during regenerative growth in this tissue.
urma predicts inhibition of the mtor molecule in the heart, kidney, and small intestine, while mtorc <dig> activation is predicted in both the kidney and small intestine, and undefined in the heart.
thus, while mtor involvement in organ regenerative growth is clear across organs, the relationships between pathway scores, molecule-level inferences, and urms are complex.
pathway maps for the nrf2-mediated oxidative stress response between fasted and 1dpf indicate consistent activation of this pathway in the kidney, liver, and small intestine .
in addition to predicted responses inferred from cpa , multiple observed genes in our dataset downstream of nrf <dig> are upregulated in these three organs, including thioredoxin , glutathione s-transferase mu  <dig> , and peroxiredoxin  <dig> , providing confirmatory evidence of nrf <dig> activation.
the response of this pathway in the heart is, however, less clear .
in the heart, nrf <dig> responses were predicted based on the observed fold-change values of only four genes, and predictions suggest inhibition of this pathway in the heart  although the direction  was not statistically significant .
it is also notable that we observed differences in the inferred consistency of integrated gene expression results and activation/inhibition inferences across organs : in the heart, only two inconsistencies are observed while the kidney, liver, and intestine have one, two, or four inconsistencies, respectively.
inferences from urma for the activation of nrf <dig> are highly consistent with activation inferences from cpa, including significant urm activation predicted for nfe2l <dig> in the liver and intestine and significant activation of nfe2l <dig> in kidney, liver, and small intestine .
in contrast, upstream regulators of this pathway were not predicted to be significantly activated or inhibited in the heart, inconsistent with the predictions given in the pathway figure .
expression response between  <dig> and  <dig> dpf
in comparison to expression between fasting and 1dpf, the ipa analyses conducted on genes differentially expressed between 1dpf and 4dpf across organs predicted a substantially smaller number of pathways as significantly enriched, the majority of which were predicted with ambiguous directions of activation.
this is likely due to the substantially smaller number of significantly differentially expressed genes identified in all organs between 1dpf and 4dpf, which is expected because 4dpf represents a sampling time intermediate between the peaking of organ growth and the regression of these phenotypes.
this time interval  aimed to capture the early stages of organs shifting expression towards organ atrophy and towards a reversion to the fasted state, and we expected to observe partial reversals in pathways predicted to be active between fasted and 1dpf, and perhaps additional new pathways involved in apoptosis and atrophy.
however, we found few consistent or clear patterns of interpretable pathway involvement between the 1dpf and 4dpf time points .
pathways predicted for this time interval include various pathways related to biosynthesis and stress response, such as unfolded protein response.
we also inferred inconsistent involvement of these pathways across organs, and none were predicted with a direction of activation .
only one pathway, mitotic roles of polo-like kinase, was predicted as significant and with a direction of activation between 1dpf and 4dpf, and was predicted only in the small intestine.
while we did infer a single lipid signaling pathway that also was indicated by cpa predictions from the fasted to 1dpf interval , the lack of predicted directions of activation and unclear involvement across organs prevents informative interpretation of the activity of this pathway between 1dpf and 4dpf.
collectively, these results suggest that the 4dpf time point may not be sufficient to capture shifts in gene expression that elucidate the mechanisms involved in the early stages of regression of organ phenotypes.
discussion
a detailed understanding of the molecular mechanisms capable of driving regenerative growth in vertebrates may provide important insights into the treatment of diverse human diseases.
because traditional vertebrate model systems offer limited insight into natural organ regenerative processes, non-traditional model systems, including snakes in general and burmese pythons in particular, hold great potential for providing unique insights into vertebrate regenerative organ growth processes.
in this study we have found that multiple integrated growth pathways, in addition to multiple stress-response pathways, appear to underlie the coordinated organ regenerative process in burmese pythons upon feeding.
despite distinct patterns of gene expression associated with growth for each organ, pathway and upstream regulatory molecule analyses reveal substantial similarities in pathways associated with post-feeding, extreme-growth responses across multiple organs.
specifically, we found evidence for a consistent interactive role of three major types of pathways underlying growth responses in python organs following feeding, including the related growth pathways mtor and pi3k/akt, lipid-signaling pathways such as ppar and lxr/rxr, and stress-response/cell-protective pathways including nrf <dig>
mtor and other growth pathways underlying organ growth
across the four organs examined, we found evidence for the involvement of the mtor signaling pathway as a key integrator of growth signals underlying post-feeding regenerative organ growth.
this pathway integrates processes for the use of energy and nutrients to regulate growth and homeostasis  <cit> .
mtor interacts with multiple other pathways, including pi3k/akt, several lipid metabolism and signaling pathways  <cit> , and the nrf2-mediated oxidative stress response  <cit>  – all of which are also active in multiple organs during growth .
mtor complex  <dig>  is the most well-characterized of the two mtor complexes and integrates signaling from growth factors, energy status, oxygen, and amino acids to promote cell growth when activated  <cit> .
the tsc1/ <dig> complex transmits upstream signals from growth factor and insulin signaling to modulate the activity of mtorc <dig> and its interaction with other pathways including pi3k/akt  <cit> .
the effector kinases of these external pathways inactivate the tsc complex through phosphorylation, thus, indirectly activating mtorc <dig>  <cit> .
akt can also directly activate mtorc <dig> through phosphorylation of an mtorc <dig> inhibitor.
in a low energy state, ampk inhibits mtorc <dig> by phosphorylating regulatory associated protein of mtorc <dig>   <cit> .
mtorc <dig> signaling is less well-understood, but is known to respond to growth factors through pi3k signaling  <cit> .
cpa of gene expression in the first 24 h after feeding indicate that involvement of the mtor signaling pathway is significant in the small intestine , but insignificant in both the heart  and kidney .
the liver lacked evidence of involvement of the mtor signaling pathway from cpa .
in urm analysis, the mtor molecule itself was predicted to be downregulated in the heart, liver, and intestine with no presence in the kidney, which contrasts our cpa results .
however, urma-predicted activation of the mtorc <dig> complex is supported in both the kidney and small intestine with undefined involvement in the heart, and the liver shows no signal for mtorc <dig> .
interestingly, cpa indicate mtorc <dig> is downregulated in the small intestine at 0–1dpf , yet this downregulated state of mtorc <dig> is based only on the downregulation of a single gene, g protein subunit beta  <dig> like , which ipa identifies as a subunit of the mtorc <dig> complex.
in contrast, ampk signaling is predicted to be downregulated in the kidney and small intestine, indicative of elevated atp levels and active mtorc <dig>  <cit>  .
it is notable that nearly all genes in the mtor pathway were associated with python orthologs that were observed as expressed across our dataset , which suggests that our inferences of non-responsive genes within the mtor pathway are biologically meaningful , rather than representative of a lack of data.
thus, mtor signaling in python tissues during regenerative organ growth may include non-canonical features compared to typical models of mtor signaling that account for the partial responsiveness of genes and targets inferred from our cpa.
our results identify mtor as a central regulator and integrator of a number of diverse growth signals that drive post-feeding regenerative organ growth in burmese pythons.
insulin signaling represents a key-regulating factor of the mtor pathway  <cit> , and we found multiple lines of evidence indicating roles of insulin signaling in post-feeding growth responses.
specifically, 0–1dpf urma inferred the activation of insr and insulin, and the inhibition of insig <dig> and insig <dig>  in the kidney, small intestine, and liver, and the inverse of these activation patterns in the heart.
insig <dig> and insig <dig> are negative regulators of scap  <cit> , which in turn regulates srebp activity.
consistent with inferences of inhibition of insig1– <dig>  urma predicted the upregulation of srebf <dig> and srebf <dig>  which provide evidence of an increase in sterol-regulatory element activity coincident with organ growth  <cit>  .
in addition to the interaction of insulin signaling and mtor activity, we also found multiple lines of evidence for pi3k/akt signaling that would interact with mtor.
our urma indicates significant downregulation of pten, an upstream regulator of the pi3k/akt pathway, across all four organs, and cpa predicts activation of the pi3k/akt signaling pathway in the small intestine and liver.
evidence from previous studies also support the role of mtor, pi3k/akt, and ampk signaling mechanisms in python post-feeding growth, at least in the heart.
western blots of python cardiac tissue post-feeding support the inference of early activation of mtor and pi3k/akt pathways by demonstrating that phosphorylated akt and mtor proteins increase significantly in abundance between  <dig> and 24 h post-feeding  <cit> .
these western blots also demonstrated phosphorylated ampk protein was upregulated within 24 h post-feeding, but lagging temporally behind the peak in phosphorylated mtor and akt  <cit> , consistent with the antagonistic relationship between ampk and mtor/akt  <cit> .
these independent lines of evidence for the roles of mtor, pi3k/akt, and ampk signaling in python post-feeding organ growth confirm our inferences of the central roles of these pathways, and support the power of pathway and urm inferences for inferring signaling mechanisms.
mapk and related pathways also appear to be prominently involved in organ growth responses post-feeding, which is sensible given their known interactions with multiple growth pathways, including pi3k/akt signaling and mtor .
our data reveal the involvement of mapk signaling most clearly in the heart, with significant enrichment and predicted inhibition of p <dig> mapk signaling and significant activation of erk <dig> signaling .
erk <dig> is a member of the mitogen-activated protein kinases  that is crucial to cell proliferation and activated in response to growth factors and oxidative stress  <cit> .
myc is a downstream transcription factor regulated by the mapk pathway and erk <dig> specifically  <cit> , and an essential regulator of development and cell proliferation .
our urma predict significant activation of myc in all four organs, indicating a broad role of active mapk signaling in post-feeding organ growth in the python.
nrf <dig> – protective function and interaction with growth pathways
one of the strongest and most consistent signals in the canonical pathway and upstream regulatory molecule analyses was the involvement of the nrf2-mediated oxidative stress response pathway.
commonly associated with anti-aging and longevity , injury repair, and mitigation of inflammation  <cit> , evidence for the central involvement of the nrf2-mediated oxidative stress response pathway in the small intestine, liver, and kidney begs the question of whether there is an important yet largely unappreciated role for stress-response signaling pathways in growth responses, and regenerative organ growth in particular.
the nrf <dig> pathway was significantly upregulated in small intestine, kidney, and liver within the first day following feeding , and the nrf <dig> transcription factor  was one of the most significant and highest in magnitude urms predicted in these three organs  .
the 24 h period following feeding in burmese pythons involves unparalleled rates and magnitudes of organ growth, and also includes massive upregulation of metabolism – up to 44-fold increases in aerobic metabolism, which is among the highest fluctuation known for any vertebrate  <cit> .
it is, therefore, sensible that activation of nrf <dig> is related to these major shifts in oxidative metabolism, and associated generation of reactive oxygen species  <cit> .
an open question, however, is what broader role the activation of nrf <dig> may play in facilitating the extraordinary growth responses associated with feeding in pythons.
for example, post-fed burmese python blood plasma has been shown to convey resistance to apoptosis to mammalian cells, even with exposure to high fatty acid concentrations that would otherwise cause cell death  <cit> ; such cell-protective qualities may be related to signals that activate nrf <dig> and/or other stress-response pathways.
interestingly, in addition to cell-protective roles of nrf <dig>  this pathway also contains multiple points of integration with various growth pathways, including those activated in python organ regenerative growth.
the nrf2-mediated oxidative stress response pathway interacts with multiple pathways predicted in our canonical pathway analysis  .
the pi3k/akt signaling pathway, predicted to be upregulated upon feeding in both the liver and small intestine, is essential for regulating the antioxidant functions of nrf <dig>  and studies have shown that inhibition of this signaling pathway leads to attenuation of nrf <dig> activities  <cit> .
this interaction is evident when examining the role of nrf <dig> in the proliferation of cancer cells.
studies have shown that nrf <dig> is able to redirect glucose and glutamine into anabolic pathways through activation of pi3k/akt signaling  <cit> .
the activated pi3k/akt pathway leads to greater accumulation of nrf <dig> in the nucleus, which allows nrf <dig> to enhance metabolic activities as well as promote cell proliferation and cytoprotection  <cit> .
the pi3k/akt signaling pathway activates mtor activity in response to growth factors, and this and previous studies  <cit>  have shown that pi3k/akt and mtor signaling are key growth pathways underlying organ regenerative growth in the burmese python.
therefore, there appears to be strong and coordinated links between growth signaling  and stress response signaling via nrf <dig> underlying organ growth in pythons following feeding.
like mtor, a large majority of genes in the nrf <dig> pathway were associated with python orthologs and were observed as expressed across our dataset , which indicates that our inferences of non-responsive genes within the nrf <dig> pathway are likely true negatives, rather than artifacts due to a lack of ortholog identification in the python.
accordingly, predicted but unobserved expression responses in the nrf <dig> pathway in pythons suggest that the absence of expected responses may represent novel or non-canonical aspects of python biology or of the organ regeneration response in pythons.
in addition to nrf2-mediated oxidative stress response, evidence for the involvement of other stress response signaling mechanisms in python post-feeding organ growth was also observed.
eif <dig> signaling, important in translational control and responsiveness to conditions of environmental stress  <cit> , is strongly downregulated in the intestine, yet, absent in the other three organs .
acute phase response signaling, which is involved in restoring homeostasis following inflammation or injury  <cit> , is predicted to be strongly downregulated in the liver and moderately upregulated .
the precise roles of these additional stress response mechanisms in regenerative organ growth in the python remains an open question, although there is strong and consistent signal for the involvement of multiple stress response pathways overall in python post-feeding organ growth.
role of lipid signaling in driving growth
previous studies have shown evidence that molecules responsible for triggering python post-feeding organ growth circulate in the blood of the burmese python  <cit> .
riquelme et al.
demonstrated that post-feeding python plasma was capable of inducing cardiomyocyte growth in pythons and mice, and that fasted python plasma supplemented with three particular fatty acids successfully stimulated cardiomyocyte growth in mice  <cit> .
because these fatty acids only facilitated a growth response in the presence of fasted burmese python serum, it is likely that python plasma contains additional factors required for successful post-feeding regenerative growth and that fatty acids are only partially responsible for stimulating growth responses.
in the heart, we found significant enrichment and predicted activation for the lxr/rxr activation pathway as well as predicted activation of this pathway  in the small intestine .
lxr is a potent activator of the srebp-1c gene  <cit> , and our data predict clear and significant activation of both srebf <dig> and srebf <dig> upon feeding in the kidney, liver, and small intestine with significant downregulation and undefined direction for srebf <dig> and srebf <dig> in the heart, respectively .
when activated, these proteins directly enhance genes important for the uptake and synthesis of various lipids.
scap, important for the activation of these sreb molecules, is also predicted to be strongly activated in the kidney, liver, and small intestine   <cit> .
we also examined ppar signaling as a potential pathway for lipid signaling during this regenerative growth, given the central role of ppar in mediating fatty acid signaling as well as its effects on gene expression  <cit> .
ppar has also been identified as an important regulator of cell survival during wound repair and regeneration  <cit> .
although cpa did not detect significant ppar signaling activation, urma significantly predicted ppara, pparg, ppargc1a, and ppargc1b involvement across organs, typically inhibited in the heart and activated in the other three organs in 0–1dpf comparisons .
given the variations in pathway and urm inferences between the heart and the other three organs, the question of whether fatty acids also play a similar stimulatory role in regenerative growth in the small intestine, liver, and kidney as they do in the heart remains.
our results do, however, argue for a poorly understood yet central role of lipid-signaling in these growth responses, and suggest that the unusually strong bioactivity of fatty acids may elicit growth through conserved canonical pathway signaling mechanisms.
early phases of organ regression following digestion
physiological studies have shown that python post-feeding organ growth peaks between 1dpf and 3dpf  <cit>  and that phenotypes begin to decline by 4dpf  <cit> .
thus, as post-feeding growth phenotypes reverse from 1dpf to 4dpf, we expected to observe shifts towards the fasted state, such as the reversal or inhibition of growth-associated pathways.
relative to comparisons between fasting and 1dpf, comparisons between 1dpf and 4dpf yielded nearly an order of magnitude fewer significantly differentially expressed genes .
accordingly, expression heatmaps  and expression profile summaries  show that expression profiles of many genes at 4dpf tend to remain elevated , or exist at intermediate levels .
we did not observe any particularly informative trends in canonical pathways and upstream regulator molecule predictions  associated with shifts in gene expression from 1dpf to 4dpf, and this result is not surprising given the relatively small number of genes that significantly change between these time points.
among the predicted pathways were several that are related to stress response and biosynthesis , although a lack of predicted direction of activation prevents detailed interpretation of the involvement of nearly all pathways predicted between 1dpf and 4dpf.
the only pathway predicted as significant and with a direction of activation between 1dpf and 4dpf was the mitotic roles of polo-like kinase pathway, which was activated in the small intestine .
it therefore remains an open question whether atrophy and other processes involved in reverting to the fasting state are controlled actively , passively  that stimulates the initial cascade of responses fades or stops), or some combination of the two mechanisms.
collectively, our results suggest that comparisons between the 1dpf to 4dpf time points may not be sufficient to predict the physiological mechanisms involved in phenotypic regression with adequate power.
further experiments, possibly with multiple later-stage time point sampling, may be required to address outstanding questions about how these growth phenotypes are reversed.
comparison of python organ regeneration to other regenerative model systems
organ regeneration in snakes represents an extreme and unique phenotype among vertebrates.
however, other examples of regenerative growth do exist among vertebrates, such as limb regeneration in salamanders  <cit> , fin regeneration in fish  <cit> , and regenerative heart growth in zebrafish  <cit>  and prenatal mammals  <cit> .
this begs the question of whether or not these regenerative responses share common mechanisms, and as we continue to better understand the mechanisms driving regenerative growth in snakes, such key comparisons can begin to be made.
while none of these other vertebrate regenerative growth systems directly parallel regenerative organ growth in snakes, regeneration of heart tissue in zebrafish is the most analogous comparison, as it occurs in adult organisms and represents regenerative growth of organ tissue specifically.
following injury or amputation of cardiac tissue, zebrafish hearts grow primarily by dedifferentiation and subsequent proliferation of cardiomyocytes  <cit> .
conversely, python hearts grow only by hypertrophy  <cit> , and therefore may be driven by largely different regenerative mechanisms.
the python small intestine, liver, and kidney, however, do grow via by hypertrophy and hyperplasia  <cit> ; while they represent different organ systems than the zebrafish heart, they may be driven by similar pathways that regulate cell proliferation in general.
indeed, there are parallels between zebrafish and python responses in the shared involvement of p <dig> mapk signaling, a negative regulator of cardiomyocyte proliferation in zebrafish  <cit>  that we infer to be inhibited in the burmese python heart between fasting and 1dpf .
additionally the mitotic roles of polo-like kinase pathway, which was the only pathway we predicted as significant and with a direction of activation between 1dpf and 4dpf  is also involved in zebrafish regenerative heart growth.
cell-cycle regulation by polo-like kinase  <dig> is an important component of cardiomyocyte proliferation in zebrafish  <cit> , and therefore may be playing a similar role in the python small intestine, although it is notable that it was not predicted as significant between fasting and 1dpf, when growth is presumably greatest in this organ  <cit> .
other pathways involved in zebrafish regenerative growth, such as igf signaling, fgf signaling, hippo signaling, and tgf-beta signaling  <cit> , were not inferred as significant based on canonical pathway analyses of either post-feeding time interval in our study of the burmese python.
tgfb <dig> and igf <dig> growth factors were, however, inferred in our urma analysis of the fasting to 1dpf interval , suggesting that there may still be some involvement of these growth factors in the regulation of regenerative growth in the burmese python.
a key conclusion based on our study is that, to our knowledge, mtor signaling and nrf2-mediated oxidative stress response pathways have not been implicated in zebrafish regenerative growth.
thus, regenerative organ growth in the burmese python appears to remain quite unique among vertebrates, both in the nature of the phenotype, and now in the molecular mechanisms underlying growth.
CONCLUSIONS
multiple coordinated growth pathways appear to play an important role in facilitating regenerative organ growth in multiple tissues of the burmese python, and the overlap of pathways across organs suggests common signaling molecules may drive this response – consistent with evidence that common factors circulating in the plasma of pythons are capable of eliciting growth  <cit> .
our analyses provide strong evidence for the involvement of particular growth and stress response pathways in post-feeding organ growth responses in multiple organs, although it is notable that our inferences of the activation versus inhibition of mechanisms was not always consistent across analyses .
as discussed above, such conflicting inferences could be due to the fundamental differences in cpa and urma , in that they are integrating very different sources of evidence, coupled with the possibility that the continuous nature of this response may survey various mechanisms during an inflection point of activity that can confound inferences of directionality.
however, contradictory inferences of mechanistic activation may also suggest that some of these core signaling pathways function differentially in snakes, or that some molecules or pathways are signaling via non-canonical mechanisms.
experiments have demonstrated that exposure to burmese python 2dpf blood serum elicits significant growth of rat cardiomyocytes  <cit> , as well as increases in size and insulin production of human pancreatic beta cells  <cit> .
these findings suggest that even if regenerative organ growth in snakes is achieved in part by non-canonical pathway or regulator activity, core aspects of signaling underlying organ growth in pythons is conserved across vertebrates.
among the most intriguing results of this study is the consistent predicted activation of the nrf2-mediated oxidative stress response pathway, and nrf2-related signaling molecules, during regenerative organ growth.
the integration of nrf <dig> signaling with other growth pathways, including mtor, provide an exciting and novel mechanistic hypothesis for how nrf <dig> and other stress-response pathways may play an important yet largely unappreciated role in regenerative growth responses in vertebrates.
additional file
additional file 1: figure s <dig>  stem analysis of all genes differentially expressed across all time points .
all significant expression profiles are shown with p-value and number of genes following that profile.
figure s <dig>  heat maps depicting activation z-scores for classes of upstream regulator molecules.
green indicates predicted activation, blue indicates predicted inhibition, white indicates that the regulator is not predicted to function in that organ, and grey indicates that the upstream regulator is predicted to have significant involvement but the activation state cannot be determined based on the gene expression data.
regulators shown on the heat maps were filtered by activation z-scores greater than | <dig> | in at least one tissue.
figure s <dig>  combined gene expression and predicted activation information for the mtor pathway in the kidney.
differentially expressed genes identified in our rnaseq data are highlighted in red  and blue  while predicted activation states are highlighted in orange  and green .
figure s <dig>  pathway prediction for the nrf2-mediated oxidative stress response in the heart.
predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set.
figure s <dig>  pathway prediction for the nrf2-mediated oxidative stress response in the kidney.
predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set.
figure s <dig>  pathway prediction for the nrf2-mediated oxidative stress response in the liver.
predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set.
figure s <dig>  pathway analysis of gene expression responses from 1dpf to 4dpf.
bar graph showing significant canonical pathways  enriched for genes differentially expressed at these time points.
pathways were filtered to include those with at least one significant p-value in one of the four organs.
bars are colored based on the predicted activation z-score for that pathway.
table s <dig>  sequencing information for all included python samples.
pe <dig> and pe <dig> stand for the sequence read type .
the year provided represents the year in which the sample was sequenced.
table s <dig>  the number of genes involved in each pathway as defined by ipa, the number of genes in the pathway that were assigned python orthologs via tblastx, and the number of those python orthologs observed with a non-zero level of expression in our dataset.
abbreviations
cpacanonical pathway analysis
dpfdays post-feeding
ipaingenuity pathway analysis
urmaupstream regulatory analysis
we thank j. castoe, r. wostl, and d. polanco for assistance with laboratory aspects of this work and the minnesota supercomputing institute and the texas advanced computing center for access to computational resources.
funding
support was provided from startup funds from the university of texas at arlington to tac, and grants from the university of texas at arlington phi sigma society to ala. additional support was provided by the national science foundation  to sms and an nih director's transformative research award  to ac.
availability of data and materials
new sequencing data as well as previously generated data are accessioned at ncbi: srp <dig>  all supporting data are included as additional files.
authors’ contributions
tac, sms, and ala conceived and designed this study.
sms performed feeding experiments and tissue sampling.
ala, tac, dcc, and drs extracted rna from tissue samples.
ala and bwp performed all gene expression, pathway, and regulatory molecule analyses.
ala, bwp and tac wrote the manuscript.
rpr, sem, and ac provided insight into physiological interpretations of analyses and contributed to manuscript sections.
all authors read, edited, and approved the final manuscript.
competing interests
the authors declare that they have no competing interests.
consent for publication
not applicable.
ethics approval and consent to participate
all animal care and tissue sampling was conducted using protocols approved by the university of alabama institutional animal care and use committee .
publisher’s note
springer nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

 the ability to massively downregulate metabolic and physiological functions during extended periods of fasting has evolved in multiple species of snakes. this downregulation of physiological form includes the atrophy of organs such as the heart, kidney, liver, and small intestine. upon feeding, the size and function of these organs, along with oxidative metabolism, is massively upregulated to accommodate digestion. of the snake species that experience these large fluctuations in physiology, the burmese python is the most well-studied. within 48 h of feeding, burmese pythons can undergo up to a 44-fold increase in metabolic rate and >100-fold increases in plasma triglyceride content. major organs also experience dramatic shifts in physiological form, including 40–100% increases in the mass of the heart, liver, pancreas, kidneys, and small intestine. this extreme organ regenerative growth and atrophy is unparalleled across vertebrates, and studies indicate that this organ growth is driven by multiple cellular processes, including cellular hypertrophy in the heart and mixtures of hypertrophy and hyperplasia in the kidney, liver, and small intestine. organ growth peaks around 1– days post-feeding, and by 10–14dpf, organ form and function, as well as gene expression patterns, have completely reversed back to fasted levels. previous studies have examined aspects of this post-feeding response using morphological and physiological assays, analyses of gene expression, and combinations of the two. together, these studies have demonstrated that transcriptional responses following feeding are extremely rapid and massive, both in the magnitude of expression changes and in the number of genes with significant differential expression. genes important in a number of developmental, metabolic, proliferative, apoptotic, and growth processes have been shown to be involved in these major shifts in organ form and function. previous studies have shown that mammalian cells respond to the growth signals in post-fed python serum, which likely indicates a conserved response to core signaling molecules. we therefore hypothesize that a relatively small number of core molecular regulatory molecules and signaling pathways may underlie these responses. however, the identification of a core set of upstream regulatory molecules and mechanisms has been hindered by the large number of genes that are significantly differentially expressed during this response, making manual interpretation of this gene expression data difficult. additionally, the lack of comparable replicated sampling across multiple organs has further prevented meaningful across-organ comparisons of changes in gene expression in previous studies. accordingly, major gaps remain in our understanding of the specific mechanisms and growth pathways that are responsible for driving these extreme shifts in burmese python organ size and function, as well as how these mechanisms may vary across different organ systems. our previous study of the burmese python feeding response addressed some gaps through the use of increased replicates and more frequent time point sampling for one organ, the small intestine. we identified over genes that were significantly differentially expressed during post-feeding regeneration in the small intestine with many of these genes being functionally linked to cellular processes such as wnt signaling, cell cycling, and apoptosis. this study also linked changes in gene expression with functional and phenotypic shifts by comparing rnaseq data with physiological and histological data. this detailed analysis was only conducted on the small intestine, however, and failed to address any upper-level signaling mechanisms and pathways. here, we leverage fully replicated organ-specific time courses detailing gene-level responses to infer canonical pathways and regulatory molecules driving post-feeding organ growth in the burmese python. we examined gene expression across four major organ systems – the heart, liver, kidney, and small intestine. we combined increased replicated sampling with statistical inferences of pathway activation and regulatory molecule prediction to identify the mechanistic drivers of cross-tissue, post-feeding organ regeneration. despite highly organ-specific gene expression responses associated with organ regenerative growth, we found evidence for high degrees of overlap in predicted pathways and regulatory molecules underlying these growth processes between organs. pathways predicted to be involved in regulating this physiological response include lxr/rxr activation, pi3k/akt, and mtor signaling. interestingly, we also found strong and consistent evidence for the involvement of nrf2-mediated oxidative stress response and other stress-response pathways in this extreme example of rapid organ growth. our results suggest that post-feeding, regenerative organ growth in the burmese python may stem from small numbers of key effector molecules mediating a core set of growth and stress-response pathways, which in turn activate diverse, tissue-specific signaling cascades. methods feeding experiments burmese pythons were obtained from commercial breeders. all animal care and tissue sampling was conducted using protocols approved by the university of alabama institutional animal care and use committee. burmese pythons were sampled at three physiological states: fasted, 1 day post-feeding and 4dpf, with the meal consumed equaling at least 25% of their body mass. previous studies have shown that organ masses and functional phenotypes climax between and dpf and that phenotypes begin to decline by 4dpf. we therefore chose sampling time points here to capture gene expression patters during the period before phenotypes climax and early in their regression. snakes were humanely euthanized by severing the spinal cord immediately behind the head, and organs were immediately extracted, snap frozen in liquid nitrogen, and stored at -80 °c. between three and six biological replicates were sampled for each time point. see additional file supplementary methods for additional details. transcriptome library generation total rna was extracted from ~50 mg of snap-frozen tissue using trizol reagent, followed by mechanical cell disruption using a tissuelyzer for 10 min at strokes/min, and precipitation of rna using isopropanol. individual illumina mrnaseq libraries were constructed using either the illumina truseq rnaseq kit or the neb next rnaseq kit, both of which included poly-a selection, rna fragmentation, cdna synthesis, and indexed illumina adapter ligation. completed rnaseq libraries were quantified on a bioanalyzer, pooled in equal molar ratios in various multiplex arrangements, and sequenced on either an illumina gaiix or illumina hiseq. quantifying and visualizing gene expression raw demultiplexed illumina rnaseq reads were quality filtered and trimmed with trimmomatic v.. in instances where the same library was sequenced in multiple different runs, reads were combined and mapped for each individual and time point. mapping of reads to the reference transcriptome of the burmese python was conducted using bwa v.. with the following parameters: mismatch penalty =  gap open penalty =  and alignment score minimum =  expression was determined using samtools v.. by counting the number of unique gene reads that mapped to an annotated transcript, while excluding reads that mapped to multiple positions. new rnaseq data for various time points and replicates was analyzed together with previously published data from other individuals and replicates. newly-generated sequencing data were archived on the ncbi short read archive. raw expression counts were normalized using tmm normalization in edger and all statistical analyses of gene expression were conducted using normalized data. we identified genes that were significantly differentially expressed between time points using two approaches. first, we estimated significant changes in gene expression between pairs of time points using pairwise exact tests for the binomial distribution calculated in edger, integrating both common and tagwise dispersion. second, to accommodate the time-series nature of the experimental design, we also conducted step-wise regression analysis of gene expression in masigpro. regression analysis enabled the detection of genes with significant patterns of differential expression across all three time points. gene expression heatmaps were generated in r and clustered with the package vegan, with gene clustering calculated using average linkage hierarchical clustering based on a bray-curtis dissimilarity matrix. we used the program stem to identify and visualize significant expression profiles for all genes in our rnaseq data. assigning homology for functional analyses to facilitate the use of various pathway activation and regulatory molecule predictions, we annotated the full burmese python transcript set with orthologous human gene ensembl identifiers. reciprocal tblastx was first conducted between anolis carolinensis and burmese python, and anolis gene ids identified as orthologous to python genes were converted to human ensembl identifiers using homology tables from ensembl’s biomart. the same process of reciprocal best blast using tblastx was performed between burmese python and gallus gallus, followed by conversion of chicken ensembl identifiers to human ensembl identifiers using homology tables from ensembl’s biomart. we also performed reciprocal best blast of the python with homo sapiens. finally, we used one-way tblastx with anolis, chicken, and human to annotate python genes that were not assigned an ortholog from reciprocal best blast. using this annotation approach, we were able to assign human ensembl ids to of total python reference transcripts. pathway and upstream regulatory molecule analyses to infer the involvement of upstream regulatory molecules and pathways, we performed core analysis in ingenuity pathway analysis, using default parameters. ipa uses gene identifiers and the fold-change value for each differentially expressed gene to identify enrichment patterns for canonical pathway analysis and upstream regulatory molecule analysis, and to infer the activation direction between particular time points. these two analyses both use observed gene expression data to infer unobserved features, but differ fundamentally in how they use expression data to make inferences. cpa predicts the involvement and activation/inhibition of canonical pathways based on observed evidence from gene expression data, specifically for genes that participate as higher-level regulatory molecules within a given pathway; analysis of observed gene expression data incorporates information from the ingenuity knowledge base to provide both a statistical value of enrichment and a prediction of the biological involvement for the pathway as a whole. in contrast, urma uses observed changes in gene expression specifically for genes at lower levels within pathways to predict activation or inhibition of regulatory molecules upstream of these genes. due to differences in these approaches, together these two methods provide a well-rounded set of comparable inferences for dissecting molecular mechanisms.fig. conceptual overview of differences between canonical pathway analysis and upstream regulatory molecule analysis. pairwise analyses on experimental gene expression data identify significantly upregulated and downregulated genes. significantly differentially expressed genes are then analyzed in two distinct ipa analyses canonical pathway analysis predicts pathway activation based on overlap of gene expression data with molecules within the pathway. d upstream regulatory molecule analysis predicts activation of specific regulatory molecules based on downstream molecules in our gene expression dataset for ipa analyses, we used only genes identified as significant in pairwise differential expression analyses between time intervals, and we input fold changes per gene averaged across biological replicates, along with our estimate of the orthologous human ensembl id for each gene. pathways important to cross-tissue physiological responses were isolated using the ipa cpa, with a right-tailed fisher’s exact test p-value of less than. we examined only those pathways that were significant, based on a predicted activation z-score, in at least one of the four organs for at least one of the post-feeding time points. for ipa analyses, the z-score is used to determine the statistical significance of the number of activated and inhibited predictions, and the sign of the value indicates the overall activation state. we used a p-value cutoff of for the cpa in ipa to reduce potentially spurious inferences. upstream regulators and hypotheses for global signaling molecules were identified using urma in ipa, with a fisher’s exact test overlap p-value threshold of. pathway network figures were modified manually from predicted network figures generated in ipa. for analysis of specific pathways, we also determined the number of genes involved in each pathway that were assigned python orthologs by our orthology analyses, and how many of these genes were expressed at some level in our dataset.  trends in gene expression across organs we used our expression data from all python samples to examine the degree to which different organ systems ‘turn on’ upon feeding and then experience ‘regression’ towards pre-feeding patterns of expression at 4dpf. we found that for each organ, the majority of differentially expressed genes showed immediate up- or downregulation from fasting to 1dpf. interestingly, each of the four organs examined appeared to experience regression towards fasting levels of expression by 4dpf to widely different extents, indicating that each organ may have its own unique temporal program of growth followed by atrophy. across organs, the heart appeared to shift towards regression the fastest. other organs experienced reversals of fasted to 1dpf expression shifts to varying degrees by 4dpf, ranging from the moderately paced small intestine and kidney, to the slow-paced liver. stem analysis further supported these temporal patterns of up-regulation and regression across organs.table numbers of differentially expressed genes between pre- and post-feeding time points for the four organs studied for each comparison, the numbers of up and downregulated genes were inferred using pairwise analysis with a benjamini-hochberg corrected p-value < regression analysis across time points, which tends to be conservative, identified hundreds of genes that were significantly differentially expressed across all three time points with genes in the heart, genes in the kidney, genes in the liver, and genes in the small intestine. of the total genes differentially expressed across all four organs, 21% are unique to the heart, 16% are unique to the kidney, 15% are unique to the liver, and 32% are unique to the small intestine. only a single gene was identified as significant in all four organs across all time points: coagulation factor x.fig. summary of significantly differentially expressed genes for all four organs identified via regression analysis. a venn diagram depicting the numbers of genes significantly differentially expressed across time points. darker colors indicate a large number of genes and lighter colors indicate a smaller number of genes. b heatmaps depicting all significantly differentially expressed genes across all time points in each organ. genes were significantly differentially expressed in the heart. there were genes significantly differentially expressed in the kidney. genes were significantly differentially expressed in the liver and genes showed significant differential expression in the small intestine to further dissect patterns of expression change following feeding, we conducted pairwise analyses of gene expression between time points for each organ. in the heart, pairwise analyses identified significantly differentially expressed genes between fasted and 1dpf, and genes were significantly differentially expressed between 1dpf and 4dpf. in the kidney, genes were significantly differentially expressed between the fasted state and 1dpf, while only eight genes were significantly differentially expressed from 1dpf to 4dpf. in contrast to the heart, we found many genes significantly differentially expressed between fasted and 4dpf in the kidney. in the liver, genes were differentially expressed within 1dpf, while only genes were significantly differentially expressed from 1dpf to 4dpf. with genes significantly differentially expressed between fasted and 4dpf, among all four organs, the liver was the least ‘reset’ to the fasting condition by 4dpf. finally, the small intestine showed higher levels of differential expression than the other three organs. within 1dpf, genes were significantly differentially expressed. from 1dpf to 4dpf, genes were upregulated and genes were downregulated, and genes were differentially expressed between fasted and 4dpf. genes and pathways implicated in differential gene expression in individual tissues to move beyond gene-specific responses and towards deciphering the mechanisms that may underlie growth responses across different organs, we identified pathways that were significantly activated/repressed between fasting and 1dpf. we found consistent evidence that the nrf stress-response pathway is activated in all tissues, except in the heart, where there was insufficient data to determine the direction of activation. we also found relatively consistent evidence for activation of the related growth pathways mtor and pi3k/akt across organs, although this inference was most significant in the heart and small intestine. we also inferred the involvement of the related pathways: lxr/rxr, lps/il-1-mediated inhibition of rxr function, ppar/rxr, and ppar signaling in multiple organs; the direction of stimulation of these pathways was both variable across organs and inconclusive in some organs. substantial involvement of cytoskeletal pathways, including actin cytoskeleton signaling and actin nucleation by arp-wasp complex, was also inferred across organs and positive in the kidney and small intestine, yet negative or inconclusive for the heart and liver, respectively.fig. canonical pathways predicted to be activated or inhibited from gene expression data. each pathway shown is significantly enriched for our genes with a fisher’s exact test p-value less than. pathways were shown only if they met our criteria for significance and had a predicted activation state in at least one organ. z-scores of indicate pathway predictions that lack a bias in the direction of gene regulation observed in our dataset. ppar signaling was also included in addition to pathway activation/repression patterns shared across organs, a number of pathways showed substantial organ-specific directionality of response. examples of this pattern include the growth-related ampk signaling pathway, erk signaling, and integrin signaling. lastly, a number of pathways appeared to be organ-specific, including p mapk and erk signaling in the heart and 14-3-3-mediated signaling in the small intestine. upstream regulatory molecule analysis of dpf responses our inferences of upstream regulatory molecules between the fasted and dpf time points supported many of the same molecular mechanisms underlying organ growth identified via cpa, such as stress response, growth, and lipid signaling pathways. we explored urm predictions for all classes of urms except biological drugs, chemicals, and micrornas. we found that many predicted urms were shared among organs, with shared among all four organs. predicted urms also showed substantial organ-specific patterns, with a large number of urms uniquely predicted for each organ. the heart showed the largest number of unique urm predictions, while only and unique urms were predicted in the kidney, liver, and small intestine, respectively.fig. predicted upstream regulators from ipa analysis of gene expression changes from fasted to 1dpf. a venn diagram of all upstream regulatory molecules analyzed. b heatmap of predicted activation z-scores for selected classes of upstream regulatory molecules. green indicates predicted activation, blue indicates predicted inhibition, white indicates the regulator is not predicted to function in that organ, and grey indicates that the upstream regulator is predicted to have significant involvement but the activation state cannot be determined based on the gene expression data. regulators shown in this heatmap were filtered by three conditions: 1) were present in at least three of the four organs, 2) are significantly predicted, and 3) have activation z-scores greater than | | in at least one organ. biological drug, chemical, and microrna categories were excluded from urm analyses to identify regulators with broadly relevant patterns across multiple organs, we focused on urms predicted significantly in at least three organs and with moderate to high activation z-score in at least one organ. a subset of the urms meeting these criteria is shown in fig. 4b, and the full set is shown in additional file 1: figure s many of these urm predictions coincided directly with predicted canonical pathways. nfe2l and atf key regulators within the nrf2-mediated oxidative stress response pathway, were predicted to be strongly activated in the small intestine, liver, and kidney, consistent with the canonical pathway analysis predictions of activation of the overall nrf pathway in these three organs. we also predicted involvement of nfkb and nfkbia, two key regulators within the nfkb signaling response pathway – this inflammatory response pathway is thought to be inhibited by activation of the nrf2-mediated oxidative stress response pathway. nfkb was predicted to be inhibited in the liver and heart, weakly activated in the kidney, and absent in the small intestine, while nfkbia was predicted to be inhibited in the liver, weakly activated in the heart and kidney, and again absent in the small intestine. activation of the growth pathways mtor and pi3k/akt were additionally supported by activation of predicted regulators such as mtorc and raf respectively, and the inhibition of pten. lipid signaling pathways such as lxr/rxr signaling, lps/il-1-mediated inhibition of rxr function, ppar/rxr, and ppar signaling were supported by several predicted urms such as pxr ligand, nr1h nr1i nr1i srebf srebf ppara, pparg, rxra, ppargc1a, and ppargc1b. these urms were consistently predicted as activated in the small intestine, liver, and kidney and either absent or predicted as inhibited in the heart. it is notable that while the inferences of activation directions of lipid signaling pathways across organs were largely ambiguous and sometimes inconsistent in our cpa, the associated urms display a consistent trend of predicted activation in the small intestine, liver, and kidney, and either predicted inhibition or absence in the heart. additionally, several urms, particularly for the mtor pathway, were predicted as inconsistent or even contradictory to the results of cpa or our experimental data. for example, while mtorc is predicted as significantly activated in the small intestine by urma, this molecule is downregulated in our experimental data. additionally, the mtor protein that is involved in forming both of the main complexes of the mtor signaling pathway is predicted to be strongly inhibited in the small intestine and weakly inhibited in the kidney and heart. both of these urm predictions appear to contradict the positive activation of the mtor signaling pathway inferred for the small intestine and heart as inferred from the cpa. in addition to urms involved in key predicted canonical pathways, upstream regulatory analysis predicted several other notable urms with strong activation or informative trends across organs. insulin and insr were both predicted as strongly activated regulators in the kidney, liver, and small intestine, suggesting a possible role of insulin receptor signaling in facilitating this regenerative response, which is also consistent with activation of the mtor pathway. myc, a regulator within the erk and p mapk signaling pathways, was predicted as activated in all four organs, although strongest in the liver. several regulators within the mapk signaling pathway were also predicted in urma, with atf and atf predicted as activated in the kidney, liver, and small intestine, and erk predicted as activated in the kidney and inhibited in the heart and liver. these urms suggest the involvement of the erk and mapk signaling pathways in this response, even though cpa predictions for these two pathways were not substantially strong. detailed dissection of nrf and mtor pathway responses to feeding we were particularly interested in our findings that the nrf stress response and the mtor growth pathways appear to be involved in post-feeding growth in multiple organs. to investigate these inferences further, we fully dissected evidence from our gene expression data for activation of these pathways by visualizing observed and inferred evidence for activation of these pathways in the context of ipa-generated pathway maps. specifically, we generated pathway predictions that integrate both observed shifts in gene expression from our data, and estimates of activation/inhibition of molecules downstream of these observed genes that are inferred based on canonical signaling patterns in these pathways. relevant to our power to detect pathway-wide signals of activity, we were able to associate over 70% of human genes within the mtor and nrf pathways with python orthologs that were expressed at some level in our dataset; thus, we expect that our power and degree of resolution of pathway activation for these particular pathways is quite good.fig. combined gene expression and predicted activation information for the mtor pathway in the heart and small intestine. a gene expression and predicted activity for the mtor pathway in the heart. b gene expression and predicted activity for the mtor pathway in the small intestine. differentially expressed genes identified in our rnaseq data set are highlighted in red and blue while predicted activation states are highlighted in orange and green. c cpa and urma results for pathways and upstream regulatory molecules involved in mtor signaling and other relevant growth pathways fig. ipa generated pathway prediction for the nrf2-mediated oxidative stress response in the small intestine. predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set. a gene expression and predicted activity for the nrf﻿ pathway in the small intestine. differentially expressed genes identified in our rnaseq data set are highlighted in red and blue while predicted activation states are highlighted in orange and green. b cpa and urma results for pathways and upstream regulatory molecules involved in nrf signaling and other related pathways pathways maps for mtor responsiveness between fasted and 1dpf show both common and divergent patterns of pathway activation among organs. the heart and kidney both show similar patterns of mtor activation, including the activation of both the mtorc and mtorc complexes. major differences in mtor activation between these two organs includes strong evidence for downregulation of ampk and the eif complex in the heart, yet, no direct and/or clear evidence for up- or downregulation of these complexes in the kidney. in the small intestine, the mtor pathway was inferred to be strongly downregulated, as is akt; ampk and the eif complex showed mixed signs of activation. it is also notable that different organs showed different levels of internal consistency in the integration of results with the known functionality within the mtor pathway. for example, the heart and kidney have either zero or one pathway connection in which gene expression results contradict the direction of activation of the pathway – for the kidney this disagreement occurs in the relationship between rsk and inhibition of tsc. in the small intestine, eight such disagreements occur, and most of these occur at the steps immediately above or below activation of mtorc and mtorc complexes. the liver was the only organ that contained no signal for the activation or repression of mtor pathway. it should be noted that inferences for mtor activation from cpa are at times contradictory to those identified via urma. while predictions based on the pathway maps indicate downregulation of mtor in the small intestine, the z-score suggests slight upregulation of this pathway during regenerative growth in this tissue. urma predicts inhibition of the mtor molecule in the heart, kidney, and small intestine, while mtorc activation is predicted in both the kidney and small intestine, and undefined in the heart. thus, while mtor involvement in organ regenerative growth is clear across organs, the relationships between pathway scores, molecule-level inferences, and urms are complex. pathway maps for the nrf2-mediated oxidative stress response between fasted and 1dpf indicate consistent activation of this pathway in the kidney, liver, and small intestine. in addition to predicted responses inferred from cpa, multiple observed genes in our dataset downstream of nrf are upregulated in these three organs, including thioredoxin, glutathione s-transferase mu, and peroxiredoxin, providing confirmatory evidence of nrf activation. the response of this pathway in the heart is, however, less clear. in the heart, nrf responses were predicted based on the observed fold-change values of only four genes, and predictions suggest inhibition of this pathway in the heart although the direction was not statistically significant. it is also notable that we observed differences in the inferred consistency of integrated gene expression results and activation/inhibition inferences across organs: in the heart, only two inconsistencies are observed while the kidney, liver, and intestine have one, two, or four inconsistencies, respectively. inferences from urma for the activation of nrf are highly consistent with activation inferences from cpa, including significant urm activation predicted for nfe2l in the liver and intestine and significant activation of nfe2l in kidney, liver, and small intestine. in contrast, upstream regulators of this pathway were not predicted to be significantly activated or inhibited in the heart, inconsistent with the predictions given in the pathway figure. expression response between and dpf in comparison to expression between fasting and 1dpf, the ipa analyses conducted on genes differentially expressed between 1dpf and 4dpf across organs predicted a substantially smaller number of pathways as significantly enriched, the majority of which were predicted with ambiguous directions of activation. this is likely due to the substantially smaller number of significantly differentially expressed genes identified in all organs between 1dpf and 4dpf, which is expected because 4dpf represents a sampling time intermediate between the peaking of organ growth and the regression of these phenotypes. this time interval aimed to capture the early stages of organs shifting expression towards organ atrophy and towards a reversion to the fasted state, and we expected to observe partial reversals in pathways predicted to be active between fasted and 1dpf, and perhaps additional new pathways involved in apoptosis and atrophy. however, we found few consistent or clear patterns of interpretable pathway involvement between the 1dpf and 4dpf time points. pathways predicted for this time interval include various pathways related to biosynthesis and stress response, such as unfolded protein response. we also inferred inconsistent involvement of these pathways across organs, and none were predicted with a direction of activation. only one pathway, mitotic roles of polo-like kinase, was predicted as significant and with a direction of activation between 1dpf and 4dpf, and was predicted only in the small intestine. while we did infer a single lipid signaling pathway that also was indicated by cpa predictions from the fasted to 1dpf interval, the lack of predicted directions of activation and unclear involvement across organs prevents informative interpretation of the activity of this pathway between 1dpf and 4dpf. collectively, these results suggest that the 4dpf time point may not be sufficient to capture shifts in gene expression that elucidate the mechanisms involved in the early stages of regression of organ phenotypes. discussion a detailed understanding of the molecular mechanisms capable of driving regenerative growth in vertebrates may provide important insights into the treatment of diverse human diseases. because traditional vertebrate model systems offer limited insight into natural organ regenerative processes, non-traditional model systems, including snakes in general and burmese pythons in particular, hold great potential for providing unique insights into vertebrate regenerative organ growth processes. in this study we have found that multiple integrated growth pathways, in addition to multiple stress-response pathways, appear to underlie the coordinated organ regenerative process in burmese pythons upon feeding. despite distinct patterns of gene expression associated with growth for each organ, pathway and upstream regulatory molecule analyses reveal substantial similarities in pathways associated with post-feeding, extreme-growth responses across multiple organs. specifically, we found evidence for a consistent interactive role of three major types of pathways underlying growth responses in python organs following feeding, including the related growth pathways mtor and pi3k/akt, lipid-signaling pathways such as ppar and lxr/rxr, and stress-response/cell-protective pathways including nrf mtor and other growth pathways underlying organ growth across the four organs examined, we found evidence for the involvement of the mtor signaling pathway as a key integrator of growth signals underlying post-feeding regenerative organ growth. this pathway integrates processes for the use of energy and nutrients to regulate growth and homeostasis. mtor interacts with multiple other pathways, including pi3k/akt, several lipid metabolism and signaling pathways, and the nrf2-mediated oxidative stress response – all of which are also active in multiple organs during growth. mtor complex is the most well-characterized of the two mtor complexes and integrates signaling from growth factors, energy status, oxygen, and amino acids to promote cell growth when activated. the tsc1/ complex transmits upstream signals from growth factor and insulin signaling to modulate the activity of mtorc and its interaction with other pathways including pi3k/akt. the effector kinases of these external pathways inactivate the tsc complex through phosphorylation, thus, indirectly activating mtorc. akt can also directly activate mtorc through phosphorylation of an mtorc inhibitor. in a low energy state, ampk inhibits mtorc by phosphorylating regulatory associated protein of mtorc. mtorc signaling is less well-understood, but is known to respond to growth factors through pi3k signaling. cpa of gene expression in the first 24 h after feeding indicate that involvement of the mtor signaling pathway is significant in the small intestine, but insignificant in both the heart and kidney. the liver lacked evidence of involvement of the mtor signaling pathway from cpa. in urm analysis, the mtor molecule itself was predicted to be downregulated in the heart, liver, and intestine with no presence in the kidney, which contrasts our cpa results. however, urma-predicted activation of the mtorc complex is supported in both the kidney and small intestine with undefined involvement in the heart, and the liver shows no signal for mtorc. interestingly, cpa indicate mtorc is downregulated in the small intestine at 0–1dpf, yet this downregulated state of mtorc is based only on the downregulation of a single gene, g protein subunit beta like, which ipa identifies as a subunit of the mtorc complex. in contrast, ampk signaling is predicted to be downregulated in the kidney and small intestine, indicative of elevated atp levels and active mtorc. it is notable that nearly all genes in the mtor pathway were associated with python orthologs that were observed as expressed across our dataset, which suggests that our inferences of non-responsive genes within the mtor pathway are biologically meaningful, rather than representative of a lack of data. thus, mtor signaling in python tissues during regenerative organ growth may include non-canonical features compared to typical models of mtor signaling that account for the partial responsiveness of genes and targets inferred from our cpa. our results identify mtor as a central regulator and integrator of a number of diverse growth signals that drive post-feeding regenerative organ growth in burmese pythons. insulin signaling represents a key-regulating factor of the mtor pathway, and we found multiple lines of evidence indicating roles of insulin signaling in post-feeding growth responses. specifically, 0–1dpf urma inferred the activation of insr and insulin, and the inhibition of insig and insig in the kidney, small intestine, and liver, and the inverse of these activation patterns in the heart. insig and insig are negative regulators of scap, which in turn regulates srebp activity. consistent with inferences of inhibition of insig1– urma predicted the upregulation of srebf and srebf which provide evidence of an increase in sterol-regulatory element activity coincident with organ growth. in addition to the interaction of insulin signaling and mtor activity, we also found multiple lines of evidence for pi3k/akt signaling that would interact with mtor. our urma indicates significant downregulation of pten, an upstream regulator of the pi3k/akt pathway, across all four organs, and cpa predicts activation of the pi3k/akt signaling pathway in the small intestine and liver. evidence from previous studies also support the role of mtor, pi3k/akt, and ampk signaling mechanisms in python post-feeding growth, at least in the heart. western blots of python cardiac tissue post-feeding support the inference of early activation of mtor and pi3k/akt pathways by demonstrating that phosphorylated akt and mtor proteins increase significantly in abundance between and 24 h post-feeding. these western blots also demonstrated phosphorylated ampk protein was upregulated within 24 h post-feeding, but lagging temporally behind the peak in phosphorylated mtor and akt, consistent with the antagonistic relationship between ampk and mtor/akt. these independent lines of evidence for the roles of mtor, pi3k/akt, and ampk signaling in python post-feeding organ growth confirm our inferences of the central roles of these pathways, and support the power of pathway and urm inferences for inferring signaling mechanisms. mapk and related pathways also appear to be prominently involved in organ growth responses post-feeding, which is sensible given their known interactions with multiple growth pathways, including pi3k/akt signaling and mtor. our data reveal the involvement of mapk signaling most clearly in the heart, with significant enrichment and predicted inhibition of p mapk signaling and significant activation of erk signaling. erk is a member of the mitogen-activated protein kinases that is crucial to cell proliferation and activated in response to growth factors and oxidative stress. myc is a downstream transcription factor regulated by the mapk pathway and erk specifically, and an essential regulator of development and cell proliferation. our urma predict significant activation of myc in all four organs, indicating a broad role of active mapk signaling in post-feeding organ growth in the python. nrf – protective function and interaction with growth pathways one of the strongest and most consistent signals in the canonical pathway and upstream regulatory molecule analyses was the involvement of the nrf2-mediated oxidative stress response pathway. commonly associated with anti-aging and longevity, injury repair, and mitigation of inflammation, evidence for the central involvement of the nrf2-mediated oxidative stress response pathway in the small intestine, liver, and kidney begs the question of whether there is an important yet largely unappreciated role for stress-response signaling pathways in growth responses, and regenerative organ growth in particular. the nrf pathway was significantly upregulated in small intestine, kidney, and liver within the first day following feeding, and the nrf transcription factor was one of the most significant and highest in magnitude urms predicted in these three organs. the 24 h period following feeding in burmese pythons involves unparalleled rates and magnitudes of organ growth, and also includes massive upregulation of metabolism – up to 44-fold increases in aerobic metabolism, which is among the highest fluctuation known for any vertebrate. it is, therefore, sensible that activation of nrf is related to these major shifts in oxidative metabolism, and associated generation of reactive oxygen species. an open question, however, is what broader role the activation of nrf may play in facilitating the extraordinary growth responses associated with feeding in pythons. for example, post-fed burmese python blood plasma has been shown to convey resistance to apoptosis to mammalian cells, even with exposure to high fatty acid concentrations that would otherwise cause cell death; such cell-protective qualities may be related to signals that activate nrf and/or other stress-response pathways. interestingly, in addition to cell-protective roles of nrf this pathway also contains multiple points of integration with various growth pathways, including those activated in python organ regenerative growth. the nrf2-mediated oxidative stress response pathway interacts with multiple pathways predicted in our canonical pathway analysis. the pi3k/akt signaling pathway, predicted to be upregulated upon feeding in both the liver and small intestine, is essential for regulating the antioxidant functions of nrf and studies have shown that inhibition of this signaling pathway leads to attenuation of nrf activities. this interaction is evident when examining the role of nrf in the proliferation of cancer cells. studies have shown that nrf is able to redirect glucose and glutamine into anabolic pathways through activation of pi3k/akt signaling. the activated pi3k/akt pathway leads to greater accumulation of nrf in the nucleus, which allows nrf to enhance metabolic activities as well as promote cell proliferation and cytoprotection. the pi3k/akt signaling pathway activates mtor activity in response to growth factors, and this and previous studies have shown that pi3k/akt and mtor signaling are key growth pathways underlying organ regenerative growth in the burmese python. therefore, there appears to be strong and coordinated links between growth signaling and stress response signaling via nrf underlying organ growth in pythons following feeding. like mtor, a large majority of genes in the nrf pathway were associated with python orthologs and were observed as expressed across our dataset, which indicates that our inferences of non-responsive genes within the nrf pathway are likely true negatives, rather than artifacts due to a lack of ortholog identification in the python. accordingly, predicted but unobserved expression responses in the nrf pathway in pythons suggest that the absence of expected responses may represent novel or non-canonical aspects of python biology or of the organ regeneration response in pythons. in addition to nrf2-mediated oxidative stress response, evidence for the involvement of other stress response signaling mechanisms in python post-feeding organ growth was also observed. eif signaling, important in translational control and responsiveness to conditions of environmental stress, is strongly downregulated in the intestine, yet, absent in the other three organs. acute phase response signaling, which is involved in restoring homeostasis following inflammation or injury, is predicted to be strongly downregulated in the liver and moderately upregulated. the precise roles of these additional stress response mechanisms in regenerative organ growth in the python remains an open question, although there is strong and consistent signal for the involvement of multiple stress response pathways overall in python post-feeding organ growth. role of lipid signaling in driving growth previous studies have shown evidence that molecules responsible for triggering python post-feeding organ growth circulate in the blood of the burmese python. riquelme et al. demonstrated that post-feeding python plasma was capable of inducing cardiomyocyte growth in pythons and mice, and that fasted python plasma supplemented with three particular fatty acids successfully stimulated cardiomyocyte growth in mice. because these fatty acids only facilitated a growth response in the presence of fasted burmese python serum, it is likely that python plasma contains additional factors required for successful post-feeding regenerative growth and that fatty acids are only partially responsible for stimulating growth responses. in the heart, we found significant enrichment and predicted activation for the lxr/rxr activation pathway as well as predicted activation of this pathway in the small intestine. lxr is a potent activator of the srebp-1c gene, and our data predict clear and significant activation of both srebf and srebf upon feeding in the kidney, liver, and small intestine with significant downregulation and undefined direction for srebf and srebf in the heart, respectively. when activated, these proteins directly enhance genes important for the uptake and synthesis of various lipids. scap, important for the activation of these sreb molecules, is also predicted to be strongly activated in the kidney, liver, and small intestine. we also examined ppar signaling as a potential pathway for lipid signaling during this regenerative growth, given the central role of ppar in mediating fatty acid signaling as well as its effects on gene expression. ppar has also been identified as an important regulator of cell survival during wound repair and regeneration. although cpa did not detect significant ppar signaling activation, urma significantly predicted ppara, pparg, ppargc1a, and ppargc1b involvement across organs, typically inhibited in the heart and activated in the other three organs in 0–1dpf comparisons. given the variations in pathway and urm inferences between the heart and the other three organs, the question of whether fatty acids also play a similar stimulatory role in regenerative growth in the small intestine, liver, and kidney as they do in the heart remains. our results do, however, argue for a poorly understood yet central role of lipid-signaling in these growth responses, and suggest that the unusually strong bioactivity of fatty acids may elicit growth through conserved canonical pathway signaling mechanisms. early phases of organ regression following digestion physiological studies have shown that python post-feeding organ growth peaks between 1dpf and 3dpf and that phenotypes begin to decline by 4dpf. thus, as post-feeding growth phenotypes reverse from 1dpf to 4dpf, we expected to observe shifts towards the fasted state, such as the reversal or inhibition of growth-associated pathways. relative to comparisons between fasting and 1dpf, comparisons between 1dpf and 4dpf yielded nearly an order of magnitude fewer significantly differentially expressed genes. accordingly, expression heatmaps and expression profile summaries show that expression profiles of many genes at 4dpf tend to remain elevated, or exist at intermediate levels. we did not observe any particularly informative trends in canonical pathways and upstream regulator molecule predictions associated with shifts in gene expression from 1dpf to 4dpf, and this result is not surprising given the relatively small number of genes that significantly change between these time points. among the predicted pathways were several that are related to stress response and biosynthesis, although a lack of predicted direction of activation prevents detailed interpretation of the involvement of nearly all pathways predicted between 1dpf and 4dpf. the only pathway predicted as significant and with a direction of activation between 1dpf and 4dpf was the mitotic roles of polo-like kinase pathway, which was activated in the small intestine. it therefore remains an open question whether atrophy and other processes involved in reverting to the fasting state are controlled actively, passively that stimulates the initial cascade of responses fades or stops), or some combination of the two mechanisms. collectively, our results suggest that comparisons between the 1dpf to 4dpf time points may not be sufficient to predict the physiological mechanisms involved in phenotypic regression with adequate power. further experiments, possibly with multiple later-stage time point sampling, may be required to address outstanding questions about how these growth phenotypes are reversed. comparison of python organ regeneration to other regenerative model systems organ regeneration in snakes represents an extreme and unique phenotype among vertebrates. however, other examples of regenerative growth do exist among vertebrates, such as limb regeneration in salamanders, fin regeneration in fish, and regenerative heart growth in zebrafish and prenatal mammals. this begs the question of whether or not these regenerative responses share common mechanisms, and as we continue to better understand the mechanisms driving regenerative growth in snakes, such key comparisons can begin to be made. while none of these other vertebrate regenerative growth systems directly parallel regenerative organ growth in snakes, regeneration of heart tissue in zebrafish is the most analogous comparison, as it occurs in adult organisms and represents regenerative growth of organ tissue specifically. following injury or amputation of cardiac tissue, zebrafish hearts grow primarily by dedifferentiation and subsequent proliferation of cardiomyocytes. conversely, python hearts grow only by hypertrophy, and therefore may be driven by largely different regenerative mechanisms. the python small intestine, liver, and kidney, however, do grow via by hypertrophy and hyperplasia; while they represent different organ systems than the zebrafish heart, they may be driven by similar pathways that regulate cell proliferation in general. indeed, there are parallels between zebrafish and python responses in the shared involvement of p mapk signaling, a negative regulator of cardiomyocyte proliferation in zebrafish that we infer to be inhibited in the burmese python heart between fasting and 1dpf. additionally the mitotic roles of polo-like kinase pathway, which was the only pathway we predicted as significant and with a direction of activation between 1dpf and 4dpf is also involved in zebrafish regenerative heart growth. cell-cycle regulation by polo-like kinase is an important component of cardiomyocyte proliferation in zebrafish, and therefore may be playing a similar role in the python small intestine, although it is notable that it was not predicted as significant between fasting and 1dpf, when growth is presumably greatest in this organ. other pathways involved in zebrafish regenerative growth, such as igf signaling, fgf signaling, hippo signaling, and tgf-beta signaling, were not inferred as significant based on canonical pathway analyses of either post-feeding time interval in our study of the burmese python. tgfb and igf growth factors were, however, inferred in our urma analysis of the fasting to 1dpf interval, suggesting that there may still be some involvement of these growth factors in the regulation of regenerative growth in the burmese python. a key conclusion based on our study is that, to our knowledge, mtor signaling and nrf2-mediated oxidative stress response pathways have not been implicated in zebrafish regenerative growth. thus, regenerative organ growth in the burmese python appears to remain quite unique among vertebrates, both in the nature of the phenotype, and now in the molecular mechanisms underlying growth.  multiple coordinated growth pathways appear to play an important role in facilitating regenerative organ growth in multiple tissues of the burmese python, and the overlap of pathways across organs suggests common signaling molecules may drive this response – consistent with evidence that common factors circulating in the plasma of pythons are capable of eliciting growth. our analyses provide strong evidence for the involvement of particular growth and stress response pathways in post-feeding organ growth responses in multiple organs, although it is notable that our inferences of the activation versus inhibition of mechanisms was not always consistent across analyses. as discussed above, such conflicting inferences could be due to the fundamental differences in cpa and urma, in that they are integrating very different sources of evidence, coupled with the possibility that the continuous nature of this response may survey various mechanisms during an inflection point of activity that can confound inferences of directionality. however, contradictory inferences of mechanistic activation may also suggest that some of these core signaling pathways function differentially in snakes, or that some molecules or pathways are signaling via non-canonical mechanisms. experiments have demonstrated that exposure to burmese python 2dpf blood serum elicits significant growth of rat cardiomyocytes, as well as increases in size and insulin production of human pancreatic beta cells. these findings suggest that even if regenerative organ growth in snakes is achieved in part by non-canonical pathway or regulator activity, core aspects of signaling underlying organ growth in pythons is conserved across vertebrates. among the most intriguing results of this study is the consistent predicted activation of the nrf2-mediated oxidative stress response pathway, and nrf2-related signaling molecules, during regenerative organ growth. the integration of nrf signaling with other growth pathways, including mtor, provide an exciting and novel mechanistic hypothesis for how nrf and other stress-response pathways may play an important yet largely unappreciated role in regenerative growth responses in vertebrates. additional file additional file 1: figure s stem analysis of all genes differentially expressed across all time points. all significant expression profiles are shown with p-value and number of genes following that profile. figure s heat maps depicting activation z-scores for classes of upstream regulator molecules. green indicates predicted activation, blue indicates predicted inhibition, white indicates that the regulator is not predicted to function in that organ, and grey indicates that the upstream regulator is predicted to have significant involvement but the activation state cannot be determined based on the gene expression data. regulators shown on the heat maps were filtered by activation z-scores greater than | | in at least one tissue. figure s combined gene expression and predicted activation information for the mtor pathway in the kidney. differentially expressed genes identified in our rnaseq data are highlighted in red and blue while predicted activation states are highlighted in orange and green. figure s pathway prediction for the nrf2-mediated oxidative stress response in the heart. predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set. figure s pathway prediction for the nrf2-mediated oxidative stress response in the kidney. predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set. figure s pathway prediction for the nrf2-mediated oxidative stress response in the liver. predicted activation state of the pathway was estimated using genes identified as significantly differentially expressed from our rnaseq data set. figure s pathway analysis of gene expression responses from 1dpf to 4dpf. bar graph showing significant canonical pathways enriched for genes differentially expressed at these time points. pathways were filtered to include those with at least one significant p-value in one of the four organs. bars are colored based on the predicted activation z-score for that pathway. table s sequencing information for all included python samples. pe and pe stand for the sequence read type. the year provided represents the year in which the sample was sequenced. table s the number of genes involved in each pathway as defined by ipa, the number of genes in the pathway that were assigned python orthologs via tblastx, and the number of those python orthologs observed with a non-zero level of expression in our dataset. abbreviations cpacanonical pathway analysis dpfdays post-feeding ipaingenuity pathway analysis urmaupstream regulatory analysis we thank j. castoe, r. wostl, and d. polanco for assistance with laboratory aspects of this work and the minnesota supercomputing institute and the texas advanced computing center for access to computational resources. funding support was provided from startup funds from the university of texas at arlington to tac, and grants from the university of texas at arlington phi sigma society to ala. additional support was provided by the national science foundation to sms and an nih director's transformative research award to ac. availability of data and materials new sequencing data as well as previously generated data are accessioned at ncbi: srp all supporting data are included as additional files. authors’ contributions tac, sms, and ala conceived and designed this study. sms performed feeding experiments and tissue sampling. ala, tac, dcc, and drs extracted rna from tissue samples. ala and bwp performed all gene expression, pathway, and regulatory molecule analyses. ala, bwp and tac wrote the manuscript. rpr, sem, and ac provided insight into physiological interpretations of analyses and contributed to manuscript sections. all authors read, edited, and approved the final manuscript. competing interests the authors declare that they have no competing interests. consent for publication not applicable. ethics approval and consent to participate all animal care and tissue sampling was conducted using protocols approved by the university of alabama institutional animal care and use committee. publisher’s note springer nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
pathways such as mtor signaling, ppar/lxr/rxr signaling, and nrf2-mediated oxidative stress response are significantly differentially regulated in multiple tissues, indicative of cell growth and proliferation along with coordinated cell-protective stress responses. consistent enrichment of the nrf2-mediated oxidative stress response indicates this pathway may be particularly important in mediating cellular stress during such extreme regenerative growth. previous studies examining post-feeding organ regeneration in the burmese python have identified thousands of genes that are significantly differentially regulated during this process. these results suggest that largely divergent patterns of post-feeding gene regulation across tissues are mediated by a core set of higher-level signaling molecules. however, substantial gaps remain in our understanding of coherent mechanisms and specific growth pathways that underlie these rapid and extensive shifts in organ form and function. here we addressed these gaps by comparing gene expression in the burmese python heart, liver, kidney, and small intestine across pre- and post-feeding time points, and by conducting detailed analyses of molecular pathways and predictions of upstream regulatory molecules across these organ systems. identified enriched canonical pathways and upstream regulators indicate that while downstream transcriptional responses are fairly tissue specific, a suite of core pathways and upstream regulator molecules are shared among responsive tissues.
 previous studies examining post-feeding organ regeneration in the burmese python have identified thousands of genes that are significantly differentially regulated during this process. however, substantial gaps remain in our understanding of coherent mechanisms and specific growth pathways that underlie these rapid and extensive shifts in organ form and function. here we addressed these gaps by comparing gene expression in the burmese python heart, liver, kidney, and small intestine across pre- and post-feeding time points, and by conducting detailed analyses of molecular pathways and predictions of upstream regulatory molecules across these organ systems.  identified enriched canonical pathways and upstream regulators indicate that while downstream transcriptional responses are fairly tissue specific, a suite of core pathways and upstream regulator molecules are shared among responsive tissues. pathways such as mtor signaling, ppar/lxr/rxr signaling, and nrf2-mediated oxidative stress response are significantly differentially regulated in multiple tissues, indicative of cell growth and proliferation along with coordinated cell-protective stress responses. upstream regulatory molecule analyses identify multiple growth factors, kinase receptors, and transmembrane receptors, both within individual organs and across separate tissues. downstream transcription factors myc and srebf are induced in all tissues.  these results suggest that largely divergent patterns of post-feeding gene regulation across tissues are mediated by a core set of higher-level signaling molecules. consistent enrichment of the nrf2-mediated oxidative stress response indicates this pathway may be particularly important in mediating cellular stress during such extreme regenerative growth. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords hyperplasiahypertrophynrf2mtorregenerationrnaseqhttp://dx.doi.org/http://dx.doi.org/ /100000154division of integrative organismal systemsiob-0466139issue-copyright-statement© the author 2017 