Feeding triggers rapid, extreme changes in intestinal form and function
We sampled small intestine tissues from 3 to 5 snakes at each of five time points: fasted, 12 h post-feeding (hrpf), 1dpf, 3dpf, and 6dpf. Within 6 days following feeding, the small intestine wet mass increases by over 100%, mucosal width doubles, and enterocyte volume doubles by 6dpf (Fig. S1a-c). At the ultrastructural level, there is a 10-fold increase in intestinal epithelial microvillus length by 3dpf (Fig. 1a). Changes in intestinal form coincide with increases in intestinal performance: nutrient uptake rates for proline, leucine, and glucose increase significantly between three- and five-fold within 12 h (Fig. 1b), and intestinal enzyme activity for aminopeptidase N (APN) and maltase increases two- to four- fold (Fig. 1c). Nutrient uptake and enzyme activity capacities significantly (p < 0.05) increase five- to seven-fold by 12hrpf (Fig. S1d-e).
Large shifts in gene expression accompany the complete remodeling of intestinal architecture and massive increase in function, with 1833 genes differentially expressed between fasted and 12hrpf (Table S1). We clustered gene expression profiles over time into seven discrete patterns (Fig. 2a). The two largest gene clusters both exhibit rapid upregulation: cluster 1 (527 genes) is characterized by sustained high expression from 12hrpf to 1dpf, and cluster 2 (466 genes) by peak upregulation at 1dpf. Functional enrichment terms associated with genes in cluster 1 include ion transport, the UPR, vesicle transport, nucleoside metabolism, and carbohydrate transmembrane transport, while cluster 2 genes are enriched for phosphatase binding, RNA polymerase II DNA binding, and coenzyme binding (Figs. 2a, S2a-c). In contrast, clusters 4 (308 genes) and 6 (268 genes) genes are significantly downregulated over the first 12hrpf and 1dpf, respectively, and both recover signaling and are upregulated above fasted levels by 6dpf. Enriched functional terms for these genes include regulation of cytokinetic process, mitotic sister chromatic segregation, and centromere complex assembly (Figs. 2a, S2d-f).
To complement our physiological and transcriptomic analyses, we quantified protein and phosphoprotein expression from a subset of samples (Table S3). Our quantitative proteomic analyses identify 922 proteins with high confidence, 286 of which are also present among differentially expressed transcripts, and we clustered proteins based on expression (Fig. 2b). Cluster 1 (331 proteins) is highly expressed through fasted and 12hrpf before steadily declining in later time points, with top enrichment for functions in metabolism, proteasomes, and glycolysis and gluconeogenesis (Fig. 2c). This is consistent with RNAseq predictions of high metabolic activity and protein degradation linked to stress in the early post-feeding response (Fig. 2a). Cluster 2 (198 proteins) is enriched for terms related to protein localization and cellular transport (Fig. 2c), with two separate peaks of expression at 12hrpf and 3dpf. The 168 proteins in Cluster 3 decrease at 12hrpf but increase thereafter, and they encompass functions related to fatty acid degradation and metabolism, starch and sucrose metabolism, and PPAR signaling (Fig. 2c). Cluster 4 (150 proteins) exhibits high levels by 1dpf and is enriched for the renin-angiotensin system, TCA cycle, and multiple metabolic pathways (Fig. 2c). In contrast, Cluster 5 (75 proteins) shows the inverse pattern of low levels at 1dpf and is enriched for cell junction assembly and organization, protein polymerization, and substrate adhesion-dependent cell spreading (Fig. 2c). Collectively, quantitative proteomic results confirm key functional inferences from RNAseq data and highlight evidence for responses to high metabolic activity and stress early in the post-feeding response. Observations of high levels of cell adhesion and junction assembly proteins and rebounding increases in the quantity of transport proteins at 3dpf suggest a second wave of major shifts in cellular function linked to decreased stress and preparation for renewed cell growth and proliferation pathways during the later phase of regeneration.
Overall, both RNA and protein quantification identify an early wave (from fasted through 1dpf) of transport, stress response, and metabolic functions. Signals at 3dpf and 6dpf highlight major changes to tissue form as cell junction organization, blood vessel development, and mitosis increase. These patterns correspond to changes measured in the intestine, where we find early peaks in uptake and metabolic activity while microvillus, cell, and tissue growth do not reach their peak levels until 3dpf and 6dpf.
Stress and growth responses modulate intestinal regenerative growth
To further investigate the functional significance of differentially expressed genes, we conducted Upstream Regulatory Molecule (URM) and Canonical Pathway activity predictions using Ingenuity Pathway Analysis (IPA; Fig. 3a-c). We also inferred the downstream effects of pathway activation patterns with the molecule activity prediction in IPA. Canonical pathways with significant activation changes group into major functional categories: uptake and metabolism, stress response, cell death, cell cycle and DNA damage response, and growth and proliferation. We find significant activity changes from three major nutrient sensing and signaling pathways throughout the post-feeding time series: Insulin Receptor Signaling, IGF-1 Signaling, and mTOR Signaling (Figs. 3a, S3, S4, S5). Two transport-associated pathways exhibit activation changes: PPARα/RXRα Signaling at 12hrpf compared to fasted and Aldosterone Signaling in Epithelial Cells at 1dpf compared to fasted (Figs. 3a, S3, S6). This is also consistent with a spike in PPAR signaling inferred by protein quantification results (Fig. 2c). Other major metabolic pathways active during the feeding response include the Superpathway of Cholesterol Biosynthesis, Oxidative Phosphorylation, and Fatty Acid β-oxidation I (Fig. 3a).
Concurrent with metabolic pathway activation, we identify broad activation of multiple stress response and stress-activated pathways, such as ERK5 Signaling, IL-8 Signaling, IL-6 Signaling, the NRF2-mediated Oxidative Stress Response, Hypoxia Signaling in the Cardiovascular System, and NF-κB signaling (Fig. 3a). At 12hrpf specifically, other activated or significantly present major stress response pathways are Protein Ubiquitination, the UPR, 14-3-3 Signaling, and AMPK Signaling (Figs. S3, S7). Based on both z-score and p-value, the NRF2 stress response is the most strongly upregulated pathway from fasted to 12hrpf (Fig. S8). Molecule activity prediction downstream of UPR activation at 12hrpf suggests suppression of apoptosis (Fig. S9), while inhibition of eIF2 Signaling at 1dpf promotes apoptosis and endoplasmic reticulum (ER) stress response (Fig. S10). Other pathways linked to apoptosis are activated at 12hrpf, including TGF-β Signaling, p38 MAPK Signaling, and Death Receptor Signaling (Fig. S3), as well as a pro-apoptosis branch of p53 Signaling (Fig. S11).
Coinciding with metabolic and stress response pathways are changes in activity of growth pathways that regulate cell death, growth, division, and differentiation (Fig. 3a). PI3K/AKT Signaling significantly increases at 12hrpf compared to fasted (Fig. S3), and IPA infers patterns of molecule activation leading to both increased cell growth and survival and decreased cell cycle progression and cell death. Other major growth and proliferation pathways that upregulate after feeding include ErbB Signaling, EGF Signaling, ERK/MAPK Signaling, and Growth Hormone Signaling (Fig. 3a). For several activated growth and proliferation pathways (e.g., JAK/STAT and PDGF), patterns of activity indicate activation of a branch of the pathway not involved in cell proliferation or even downstream suppression of proliferation (Figs. S12, S13). Many cell cycle-specific pathways are inhibited at 12hrpf (Fig. S14): ATM Signaling, Role of BRCA1 in DNA Damage Response, PTEN Signaling, and Cell Cycle: G2/M DNA Damage Checkpoint Regulation. Additionally, although activation of Cell Cycle Control of Chromosomal Replication is significant but nondirectional (i.e., z-score could not be calculated), downstream molecule activity indicates inhibition of the origin recognition complex and DNA replication (Fig. S15). Downstream from these pathways, IPA predicts inhibition of G1/S phase transition, homologous recombination, DNA repair, checkpoint control, and cell cycle progression.
To investigate whether early stress and growth signaling responses cooperate to regulate tissue growth under stress, we examined the relationships and overlap among pathways with activity changes during the period from fasted to 12hrpf (Fig. 2b). We identify a large cluster of genes shared among pathways involved in growth, differentiation, and proliferation that functionally overlap and interact with the NRF2 and 14-3-3-mediated signaling pathways (Fig. 2b). These findings highlight overlapping signaling connections that mechanistically link stress response pathways (NRF2 and 14-3-3) with growth pathways (including Insulin/mTOR, JAK/Stat, ERK/MapK) that are co-stimulated within 12 h following feeding (Fig. 2b). Overall, canonical pathway analyses highlight an early increase in transport and high-level metabolic pathways followed by a major surge of stress response and cell death mechanisms. These coordinate to suppress cell cycle progression and prevent significant proliferation while the tissue is under extreme metabolic stress, and near the end of the regeneration timeline, these pathways downregulate as cell proliferation and growth resume.
Identification of essential regulatory molecules and kinases directing regeneration
Based on evidence for the interaction of key stress response and growth pathways, we investigated specific regulatory molecules that might mediate these relationships to control the regenerative response. In analysis of RNAseq data across timepoints, the top upregulated URMs are PGDF BB, a regulator of vascular permeability ; XBP1, a high-level regulator of the UPR; ATF4, a transcription factor activated by both mTORC1 and integrated stress response pathways ; and NFE2L2, a top regulator of the NRF2-Mediated Oxidative Stress Response (Fig. 3c). From fasted to 12hrpf, 78 of 95 target molecules for XBP1 are consistent with activation of XBP1, with top targets including ApoA1 and HSPA5 (BiP) (Table S2), and 100 of 137 target molecules are consistent with activation of NFE2L2 (z-score 6.448, p < 0.01) (Tables S3). These two URMs are inactivated and inhibited, respectively, by 6dpf. Several growth and proliferation URMs have the highest activity at 3dpf or 6dpf, near the end of the regenerative response: INSR, EGF, ERK, TGFB1, NF-κB, MITF, ELK1, and GH (Fig. 3c). One early activated URM is RAF1, a known oncogene [28, 29] but which is also capable of tumor suppressant activity related to decreased IL-6 and JAK/STAT3 signaling activity , and simultaneous IL-6 and STAT3 activation has been shown to promote tumor formation in gastrointestinal tissues . Multiple toll-like receptors (TLRs) – TLR3, TLR4, and TLR9 – are significantly activated at 6dpf (Fig. 3c). TLRs promote proliferation and cell growth in the intestine; for example, TLR4 expression specifically induces EGF programs of proliferation . Although NF-κB pathways are activated from 12hrpf, as a URM NF-κB is only significantly activated at 6dpf, alongside Ap-1. Both are important transcription factors stimulated by TLRs .
Our phosphoproteomic analyses identify 222 phosphoproteins with quantitative changes across time points, which we group into four clusters comprising 114, 67, 38, and 3 phosphoproteins, respectively (Fig. 3d). Cluster 1 increases abundance by 12hrpf, indicating an important role for post-translational modifications in early regenerative responses. To understand mechanistic links between phosphoprotein levels and kinases that modulate these responses, we conducted kinase enrichment analyses (Fig. 3d-e). These analyses implicate activity by the kinases CSNK1E, CSNK2A2, PRKCB, RPS6KA3, GSK3β, CSNK2A1, CDK2, and AKT1 (p < 0.05) for Cluster 1 (Fig. 3e). Gene set enrichment analyses indicate that these phosphoproteins function in mRNA splicing and processing, chaperone protein regulation, and the IRE1-mediated UPR (Fig. 3g). Cluster 2 proteins are generally dephosphorylated after feeding until 3dpf (Fig. 3d). Like cluster 1, proteins in this cluster are enriched for the kinases CSNK1E, RPS6KA3, and PRKCB, and additionally for CAMK2, PRKACB, MLYK, MAPK8, MKNK1, SGK1, and SRC (Fig. 3f), indicating loss of phosphorylative activity of these kinases. Enrichment analyses place these phosphoproteins in pathways for smooth muscle contraction, platelet activation, and apoptosis (Fig. 3h). The two smallest phosphoprotein clusters both peak at 3dpf (Fig. 3c) and were analyzed together for enrichment. These clusters are also enriched for CSNK1E in addition to ATM, MAPK14, PRKAR2B, and PKD2.
Temporal coordination of growth and stress response pathway networks
To further investigate the interconnected nature of signaling networks involved in the regenerative response and the temporal dynamics of these interactions, we visualized overlapping pathways activated across all time points (Fig. 4a). All the top centralized nodes (Freeman degree ≥60) except one occur at 1dpf, highlighting 1dpf as a prominent peak of integrated signaling (Fig. 4a). These pathways include JAK/STAT Signaling, Renin-Angiotensin Signaling, Molecular Mechanisms of Cancer, IL-8 Signaling, CXCR4 Signaling, LPS-Stimulated MAPK Signaling, EGF Signaling, PEDF Signaling, NF-κB Activation by Viruses, and B Cell Receptor Signaling (Fig. 4a). Other centralized nodes of pathway activity shift throughout the regenerative response. Although multiple key nutrient sensing and stress response pathways are active at 12hrpf, they do not have high centrality to the broader regenerative network; this suggests that they do not have direct mechanisms of feedback or interaction that link to the central signaling network hub (Fig. 4a). At 3dpf, we find that centralized nodes include tumor-relevant pathways such as UVC-Induced MAPK Signaling and HGF Signaling (Fig. 4a). Collectively, our results identify 1dpf as a major temporal hub of signal activity and integration underlying the regenerative response and highlight the extensive degree of interactions among pathways that coordinate the regenerative response, with important focus on IL-8 Signaling and CXCR4 Signaling at 1dpf and MAPK signaling pathways through 3dpf (Fig. 4b).
Integrated inference of a novel signaling hypothesis driving intestinal regeneration
Integrated evidence from physiological, RNAseq, proteomic, and phosphoproteomic analyses suggest a model of mechanistic interactions among signal regulators and pathways that modulate regeneration in the snake intestine post-feeding that we organize into three major temporal phases of activity (Fig. 4c). The first largescale change is massive nutrient uptake, triggering response in transporter activity, RAAS activation, and PPARα/RXRα signaling, matched by extreme metabolic activity increase (Fig. 4c). Our results illustrate that nutrient uptake and metabolic shifts activate a major axis of insulin/mTOR signaling and PI3K/AKT/PTEN signaling indicated by increased insulin, IGF-1, and mTOR signaling and suppression of ATM and PTEN signaling. PTEN antagonizes PI3K/AKT activity , while mTOR is a negative regulator of ATM signaling, halting cell cycle checkpoints and sensitizing cells to oxidative stress  and apoptosis  (Fig. 4c). Changes in PI3K/AKT signaling along this axis modulate cell cycle arrest and loss of stability in adherens junctions, while uptake, metabolism, and increased transcriptional activity generate extreme ER and oxidative stress responses. These cellular stress responses are presumably cell-specific and heterogenous throughout the tissue, supported by evidence for both pro-survival and apoptotic signaling regimes moderated by the UPR and 14-3-3 signaling pathways. PRKCD, which we observe at 12hrpf, directly phosphorylates NRF2 to promote antioxidant activity , and NRF2 and AMPK signaling accompany other stress responses and likely support pro-survival stress responses in cells that survive, while TGF-β, CASP3, and p38 MAPK regulate the apoptotic responses.
At 12hrpf, our results show that regulatory kinase activity includes phosphorylation by CSNK2A1 and CSNK2A2 (Fig. 3e-f), two key subunits of CK2, a pleiotropic protein kinase that functions in many signaling networks, including Wnt and the PI3K/AKT/PTEN axis . CK2 promotes survival through enhanced NF-κB activity to regulate the inflammatory response [39, 40], thus activating transcription of both IL-6  and IL-8 , which are activated in the boa intestine (Fig. 3a-b). IL-8 Signaling is a central regulator with high overlap in expressed pathways at 1dpf and 3dpf (Fig. 2d), and both IL-8 and IL-6 signaling pathways are significantly active throughout the time series (Figs. 3a and 4a). Inhibition of PKA activity, observed in decreased activity of subunits PRKACB and PRKAR2B in the early post-feeding response, is associated with increased IL-6 signaling . At this time, RAF1 activity may suppress JAK/STAT3 signaling, protecting from tumorigenesis in a high IL-6 and IL-8 signaling environment.
After the major axis of CK2-driven activity at 1dpf in IL-8, the regenerative response switches to suppress stress response pathways and allow proliferation and growth to increase. Although inhibition of eIF2 signaling in the early response is connected to increased UPR activity, the increasing inhibition of eIF2 signaling at 3dpf and 6dpf may in turn suppress the UPR . Additionally, GSK3β activity changes at 12hrpf and 1dpf are consistent with downstream inhibitory effects on NRF2 stress response . Finally, TLR signaling at 6dpf promotes late-phase NF-κB and Ap-1 expression to drive proliferation.
Conserved regenerative signaling mechanisms across divergent snake models
To test whether multiple snake lineages with regenerative capacity utilize the same underlying mechanisms, we compared our results to those from a previous study that analyzed other snake lineages that do (pythons and rattlesnakes) and do not (water snakes) regenerate following feeding . We evaluated responses to feeding across species using fasted vs. 1dpf RNAseq-based inferences of canonical pathway activation and URM activity. Across species, we find distinct cyclic temporal progressions of gene expression based on a generalized PCA, with clusters stratified by species indicative of species-specific nuances in gene expression (Fig. 5a). Among significantly differentially expressed genes between fasted and 1dpf for each species, boa and python share 285 genes in total, while over half of the genes differentially expressed in the boa are unique to that species (Fig. S16). Differences in both the temporal rate of the response  but in the number of RNAseq reads per species may contribute to some of this observed variation among species (Table S5). Genes uniquely expressed in boas and pythons function in metabolic processes, protein localization, and the UPR (Fig. 5b), and genes unique to all regenerating species are involved in COPII-vesicle transport and protein folding and modification (Fig. 5c).
Inferences of URM and canonical pathway activation highlight many signaling patterns that are shared among regenerating species with absent or opposing activation in the non-regenerating water snake. These results implicate a core set of pathways and molecules central to the onset (1dpf) of the regenerative response across species (Figs. 4b-c). Only the regenerating species show URM activation of IL24, HSF2, NFE2L1, and the PI3K family, and inhibition of MAP 4 K4 and ubiquitin, and the boa and python both demonstrated activation of LARP1 (Fig. 5d). Interestingly, the URMs XBP1 and NFE2L2 associated with UPR and NRF2 stress response pathways are highly activated in all species, while other stress response URMs are uniquely activated only in the three regenerating species (HSF2, NRF1, NRF2L) (Fig. 5d). CPA comparisons also indicate that all species significantly activate UPR and NRF2-Mediated Oxidative Stress Response, and the Superpathway of Cholesterol Biosynthesis (Fig. 5). Importantly, several pathways with shared activation patterns among regenerating species are distinct from the water snake and may represent conserved features of the snake regenerative response (Fig. 5e).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.