Thermal stress induces tissue damage and a broad shift in regenerative signaling pathways in the honey bee digestive tract

Honey bee colonies in the USA have suffered from a higher than usual rate of die-off in recent years, going from a pre-2006 expected loss rate of 15% (van Engelsdorp et al., 2008) to over 40% in 2015–2016 (Kulhanek et al., 2017). There seems to be no single cause for this increase in honey bee disease. Rather, stresses that include nutritional deficiencies due to loss of appropriate forage, chemical poisoning from pesticides, changes to normal living conditions brought about through large-scale beekeeping practices, myriad environmental changes as a result of climate change, and parasitism by arthropod pests and pathogenic microbes likely work in concert to cause tissue pathology, disease and mortality (Steinhauer et al., 2018). In seeking to understand how stresses might synergize to impact honey bee health, efforts have been undertaken to more completely define common cellular processes and cell stress pathways that are impacted by multiple stressors.

With a changing climate, an increase in the frequency of severe weather events, such as heat waves, is anticipated (Diffenbaugh et al., 2017). Forming an understanding of the effects of climate change on various organisms is of paramount importance (Tomanek, 2010). Exposure to temperatures outside of optimal ranges affects physiological processes and activities, including locomotion and reproduction, with implications for individual survival (Stillman, 2019). The ability to function normally outside of ideal ranges is likely a critical predictor of the success of a species during climate change. Thus, strategies that limit damage caused by thermal stress provide thermotolerance. At the cellular level, thermotolerance is thought to depend largely on the capacity to protect the proteome from thermal stress through the actions of pathways of the proteostatic network (Klaips et al., 2018), including the heat shock response (HSR). However, thermotolerance likely also includes the facility to limit the resulting damage occurring at the tissue level, organismal level and, in social insects, colony level.

Individual honey bees appear to have a high capacity to endure thermal stress (Free and Spencer Booth, 1962; Severson et al., 1990; Elekonich, 2008; Abou-Shaara et al., 2012; Kovac et al., 2014). The temperature of individual bees can increase significantly above steady state to levels that are dangerous to other organisms. For example, the temperature of individual forager bees can reach up to 49°C in flight (Elekonich and Roberts, 2005) and bees involved in endothermic shivering might be expected to reach high temperatures as well. One reason for this high-level endurance is likely their robust HSR (Sala et al., 2017), which we previously characterized in workers (McKinstry et al., 2017) and queens (Shih et al., 2020).

However, less is known about other mechanisms of thermotolerance operating at the tissue and organismal levels in this or other species. Seeking to further characterize the response to thermal stress in honey bees, we used RNAseq to identify novel changes in gene expression after heat stress. We chose midgut tissue as it represents a hub for interactions between infectious disease, nutrient acquisition and regulation, and environmental toxin exposure. We found a dramatic induction of signaling pathways that have been found to be associated with regeneration of the digestive tract in the fruit fly, including the Hippo and JAK/STAT pathways. We also observed cell death and shedding of epithelial cells, which likely initiates this apparent regenerative program. Alterations in some of these same pathways in other tissues after thermal stress may indicate the existence of a shared program of response following this type of insult. This study provides important foundational characterization of the damage response and regeneration program in this key pollinating species. In addition, our research suggests that a strong tissue regeneration program may contribute to thermotolerance at the tissue and organismal level. Further work is required to determine the relative importance of tissue regeneration in honey bee thermotolerance.

Honey bees were collected from the landing board of outbred colonies in New York, NY, USA, consisting of a typical mix of Apis mellifera subspecies found in North America, at different times during the months of April to October. Only visibly healthy bees were collected and all source colonies were visually inspected for symptoms of common bacterial, fungal and viral diseases of honey bees. After cold anesthesia, the following tissues were dissected for gene expression analysis: head tissue (predominantly brain, sensory organ tissue and hypopharyngeal glands), midgut, thorax tissue (predominantly flight muscle) and abdominal wall tissue (predominantly fat body). Tissues were placed in RNAlater (Invitrogen, San Diego, CA, USA) for storage prior to RNA extraction.

For all caged experiments, honey bees (8–10) were selected as above and kept in 177.4 ml (6 fl oz) Square-bottomed Drosophila Stock Bottles (VWR) plugged with modified foam tube plugs (Jaece Industries). Bees were maintained in incubators at 35°C (unless otherwise stated) in the presence of PseudoQueen (Contech, Victoria, BC, Canada) as a source of queen mandibular pheromone (QMP) and used as per the manufacturer's instructions. For heat shock, bees were maintained for 4 h in cages at 35 or 45°C, and were not kept in cages longer than 24 h. Bees were fed 33% sucrose via a modified 1.5 ml screw-cap tube.

RNA was prepared from bees from the described tissues by manually crushing the tissue of interest with a disposable pestle in Trizol reagent (Invitrogen, San Diego, CA, USA) and extracting the RNA as per the manufacturer's instructions. RNA was subsequently DNase I treated by RQ1 RNase-Free DNase (Promega, Madison, WI, USA) and quantified. cDNA was synthesized by reverse transcription using approximately 1 µg of RNA with the iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA). Typically, 1 µl of cDNA was then used as a template for quantitative PCR (qPCR) to determine the expression levels of genes of interest using the iQ SYBR Green Supermix (Bio-Rad) in an iCycler thermo-cycler (Bio-Rad). Primer sequences for qPCR developed for this study are in Table S1. Primer sequences for the reference gene β-actin were from van Englesdorp et al. (2009). The difference between the cycle threshold number for β-actin and that of the gene of interest was used to calculate the level of that gene relative to β-actin using the 2^(−ΔCT) method (Schmittgen and Livak, 2008). All qPCR data were confirmed in at least 3 independent experiments.

For RNA-Seq analysis, RNA was isolated from midgut samples as above. After quantifying and checking the purity of the RNA, it was shipped to GENEWIZ (South Plainfield, NJ, USA). RNA was subjected to additional quality control analysis before mRNA was enriched by poly-A selection for library preparation using the NEBNext Ultra RNA library preparation kit. Six individual bee RNA-Seq libraries (3 from each group) were prepared. Sequencing was then performed using the paired end 150 bp sequencing configuration on the Illumina HiSeq 4000 platform. Sequence reads were trimmed to remove possible adapter sequences and nucleotides with poor quality using Trimmomatic v.0.36 (Bolger et al., 2014). The reads were then mapped to the Apis mellifera reference genome then available on NCBI (Amel 4.5 version) using the STAR aligner (Dobin et al., 2012). The RNA-Seq aligner is executed using a splice aligner which detects splice junctions and incorporates them to help align the entire read sequences. BAM files were generated as a result of this step. Gene hit counts were calculated using Feature counts from the Subread package. After mapping and total gene hit counts calculation, the total gene hit counts table was used for downstream differential expression analysis using DESeq2. Genes with an adjusted P-value <0.05 and absolute log₂-fold change >1 were designated significant differentially expressed genes for each comparison (Dataset 1, deposited in Dryad: https://doi.org/10.5061/dryad.1ns1rn8t6). The RNA sequence information in this study has been submitted to the Gene Expression Omnibus database under accession number GSE159083 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159083). Volcano plot analysis (Fig. 1B) shows the global transcriptional change between groups. All the genes are plotted and each data point represents a gene. The log₂-fold change of each gene is represented on the x-axis and log₁₀ of the P-value on the y-axis. An adjusted P-value of 0.05 and a log₂-fold change of 1 are indicated by red dots. These represent upregulated genes. And an adjusted P-value of 0.05 and a log₂-fold change of −1 are indicated by blue dots. These represent downregulated genes. Differentially expressed genes (DEG) are reported as LOC Gene IDs found in NCBI for both the Amel 4.5 version and more current Amel_HAv3.1 version. These LOC Gene IDs also have corresponding GB numbers in HymenopteraMine. To further investigate the function of the DEGs, we used a resource kindly provided by Dr Michelle Flenniken (Brutscher et al., 2017) in which they used reciprocal BLAST+ (Camacho et al., 2009) to identify the Drosophila melanogaster orthologs and homologs of all honey bee genes (Elsik et al., 2014). Gene ontology (GO) analysis were performed with DAVID using both the honey bee gene lists (from Amel 4.5 genome version) and the lists of the corresponding fruit fly homologs (Dennis et al., 2003; Huang et al., 2009), and GO accessions were mapped onto the corresponding pathways using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa and Goto, 2000). The predicted structure of the putative unpaired (UPD) proteins was determined using Porter 4.0 and PaleAle 4.0 (http://distillf.ucd.ie/porterpaleale/quickhelp.html) and the presence of signal sequences was predicted using Signal 4.1P (http://www.cbs.dtu.dk/services/SignalP/).

For histological analysis, preparation and staining of slides was performed by Laudier Histology, New York, NY, USA. Briefly, samples were removed from stabilization liquid (hand sanitizer), rinsed well with distilled water and placed into 10% alcoholic zinc formalin fixative for 48 h. Post-fixation, samples were dehydrated with a combination of 2-ethoxyethanol and acetonitrile. Samples were then infiltrated with several changes of a custom hydrophobic acrylic resin (base reagent methyl methacrylate, Sigma-Aldrich, St Louis, MO, USA) and subsequently polymerized into tissue blocks. Polymerized blocks were subsequently sectioned on a rotary microtome at 3 μm with a tungsten carbide knife and cut sections were placed onto adhesive glass and dried. Sections were deplasticized and stained with a modified Hematoxylin and Eosin stain (using chromium mordant instead of the commonly used aluminium sulfate or potassium sulfate) and coverslipped with an acrylic-based mounting media. Slides were examined at 20× magnification using a NIKON Eclipse E600FN microscope (Nikon, Melville, NY, USA) with brightfield illumination.

For cryosectioned samples, midguts were processed with 4% paraformaldehyde (PFA) for 1 h, followed by treatment with 5% and 20% sucrose/PBS, then embedded in optimal cutting temperature (OCT) compound (Thermo-Fisher Scientific, Waltham, MA, USA) and stored at −80°C. Midguts were sectioned at 8 μm using a Leica CM1860 set to −25°C. The slides were heated at 37°C for 30 min and stored at −80°C.

Slides were washed in PBS, 4% PFA and 0.1% Triton X-100/PBS to rehydrate and permeabilize the sections. During subsequent incubation steps, slides were marked with a hydrophobic pen and placed in a humidified container to minimize evaporation; the use of coverslips was avoided to preserve tissue morphology. Slides were then treated with TdT enzyme, TdT reaction buffer and 5-ethynyl-2′-deoxyuridine (EdU) triphosphate as per the manufacturer's instructions (Life Technologies) and incubated at 37°C for 60 min. Slides were washed in 3% BSA/0.1% Triton X-100/PBS solution, then treated with Click-iT Plus TUNEL Assay reagent with Alexa 594 (Life Technologies) as per the manufacturer's instructions at 37°C for 30 min. The slides were mounted with VectaShield (Vector Laboratories, Inc.) and left to drain for 30 min before microscopy. Slides were examined at 20× magnification using a NIKON Eclipse E600FN microscope (Nikon, Melville, NY, USA) with DAPI, FITC and TEXAS RED Filter cube.

For a complete list of genes encoding UPDl proteins in bee genomes, we used each honey bee protein coding sequence as the query to search other bee genomes for which non-redundant sequences are available (Apis dorsata, Apis cerana, Apis florea, Osmia lignaria, Osmia bicornis, Eufriesea mexicana, Megachile rotundata, Habropoda laboriosa, Nomia melanderi, Melipona quadrifasciata, Ceratina calcarata, Dufourea novaeangliae, Bombus impatiens, Megalopta genalis, Bombus vancouverensis, Bombus bifarius, Bombus terrestris, Bombus vosnesenskii and Frieseomelitta varia (Table S2). We also included the D. melanogaster UPD proteins and Homo sapiens LEPTIN in further analysis. Proteins from A. cerana and D. novaeangliae were not used in the subsequent analysis because all three UPDl proteins could not be found in these genomes. Protein alignments were generated with MUSCLE (Edgar, 2004) using default parameters and inspected manually. A maximum-likelihood phylogenetic tree was inferred from the resulting alignment using the RaxML program (version 8) (Stamatakis, 2014) with the GAMMA model. Bootstrapping was conducted with 100 replicates. The resulting tree was visualized and annotated using ggtree in R.

For analysis, data were log₁₀ transformed and compared using unpaired t-tests with Welch's correction when values fitted normal distributions or Mann–Whitney U non-parametric tests when they did not fit normal distributions. Normality was assessed using Shapiro–Wilk tests. When more than two groups were being compared, data were compared using one-way ANOVA with Tukey's multiple comparison test when values fitted normal distributions or a Kruskall–Wallis test.

To use an unbiased approach to identify novel genes induced by thermal stress in the honey bee digestive tract, we performed transcriptome profiling (RNASeq) of midguts from bees maintained at either 35 or 45°C for 4 h. There was no difference in the survival of honey bees exposed to these different temperature conditions (Fig. 1A), but the higher temperature is known to induce thermal stress and activate the HSR (McKinstry et al., 2017). We chose midgut tissue as it represents a hub for interactions between infectious disease, nutrient acquisition and regulation, and environmental toxin exposure. We used bees from the landing board of the colony, which are most likely foragers, because they typically come into direct contact with more stressors than nurse bees (e.g. Vannette et al., 2015) and are therefore likely to have the most robust stress responses. Transcriptome analysis revealed that HSR induction increased the expression of 1393 genes as compared with control bees and decreased the expression of 535 genes compared with control bees (Fig. 1; Dryad Dataset 1: https://doi.org/10.5061/dryad.1ns1rn8t6), many of which are not well characterized. We performed a functional enrichment analysis on the upregulated and downregulated genes and mapped these onto KEGG pathways. Genes from the following pathways were overrepresented: upregulated genes – Endocytosis (ame04144), Protein processing in endoplasmic reticulum (ame04141), Inositol phosphate metabolism (ame00562), Hippo signaling pathway – fly (ame04391), Wnt signaling pathway (ame04310), TGF-beta signaling pathway (ame04350) and mTOR signaling pathway (ame04150) for both Apis (Table S1) and Drosophila (Table S1); downregulated genes – Non-homologous end-joining (ame03450), Mismatch repair (ame03430) and DNA replication (ame03030). Importantly, the genes we found to be upregulated by heat shock in bees in our previous study (McKinstry et al., 2017) were also found to be increased in this RNAseq experiment.

We confirmed these findings with bees maintained at either 35°C (n=7) or 45°C (n=10) for 4 h using RT-qPCR. As before, the hsp70 chaperone gene encoding Hsp70Ab was highly upregulated by heat shock (Fig. 2). Interestingly, Hsp83l was shown to be upregulated by RNAseq, differing from our previous result (McKinstry et al., 2017). Examining our previous primers, we found that a new gene model was available with the new version of the honey bee genome. Designing a new primer set, we did find that Hsp83l was induced by heat shock (Fig. 2). Relative to β-actin, we observed robust induction of two co-chaperones, Activator Of Heat Shock 90 kDa Protein ATPase Homolog 1 (Ahsa1) and BAG family molecular chaperone regulator 2 (Bag2) (Fig. 2). In addition, a number of Hippo pathway associated genes (Kibra, Myc and Wengen) and JNK pathway associated genes [Jun-related antigen (Jra), Hemipterous (Hep), Kayak (Kay) and Puckered (Puck)] were all increased after thermal stress (Fig. 2; Table S2). The relative expression of all genes tested had returned to similar levels in bees maintained at 35 or 45°C for 4 h by 20 h after cessation of the thermal stress (Fig. S1). We previously showed that β-actin levels were similar irrespective of temperature as assessed by cycle threshold values (McKinstry et al., 2017).

We quantified induction of these damage response genes in head tissue (predominantly brain, sensory organ tissue and hypopharyngeal glands, n=6), midgut (n=6), thorax tissue (predominantly flight muscle, n=6) and abdominal wall (predominantly fat body, n=6 for 35°C and n=5 for 45°C) tissues at 35 and 45°C for 4 h. Relative to β-actin, we observed robust induction of the hsp70 chaperone genes encoding Hsp70Ab and Hsc70-4, the Hsp90 chaperone encoding gene Hsp83l, and genes for the two co-chaperones, Ahsa1 and Bag2, in all tissues. The Hippo pathway associated gene Myc and JNK pathway associated gene Jra were also induced in all tissues. Other genes, including Kibra, Hep, Kay and Puck, were increased in multiple tissues after thermal stress, while Wengen was only upregulated in the digestive tract (Fig. S2). The magnitude of induction differed between tissues.

One of the most highly upregulated genes (gene ID 102655202, induced 164.5-fold), encodes a protein with a low level of similarity to Drosophila JAK/STAT ligand Unpaired 2 (Upd2). A homology search using this amino acid sequence revealed two other putative honey bee Upd genes, gene ID 100577882, which was upregulated 75.3-fold (Table S2) and gene ID 100577920, which was not altered in expression (Table S2). As with the D. melanogaster Upd genes (Hombría et al., 2005), all three putative honey bee Updl genes are located in a single cluster and we assigned the names Upd-like A–C to these genes based on their location within the gene cluster (see below). The similar organization between the two species suggested the conservation of this cluster in the last common ancestor. However, our analysis of neighboring genes did not show evidence that the gene clusters from the two species were syntenic, perhaps discounting this possibility. Using qPCR, we validated heat-dependent induction of these putative Upd-like genes after thermal stress in midgut tissue from bees maintained at 35 or 45°C for 4 h. Both gene ID 102655202 (UpdlC) and gene ID 100577882 (UpdlB) were highly upregulated by heat shock in midgut tissue, while gene ID 100577920 (UpdlA) was not (Fig. 3A,B). We examined induction of the three Updl genes in head tissue (predominantly brain, sensory organ tissue and hypopharyngeal glands, n=6), midgut (n=6), thorax tissue (predominantly flight muscle, n=6) and abdominal wall (predominantly fat body, n=6 for 35°C and n=5 for 45°C) tissues at 35 and 45°C for 4 h. Both UpdlC and UpdlB were increased in multiple tissues after thermal stress while UpdlA was only upregulated in head tissue (Fig. S1).

All three putative honey bee orthologs share some regions of marked similarity with all three Drosophila UPD proteins (Fig. S3), including a WxxxC motif, a W/FJP motif and a LC-containing motif (Fig. 4A) containing the C residue thought to be conserved between invertebrate UPD proteins and vertebrate LEPTIN proteins (Londraville et al., 2017). To explore conservation of these proteins in bees and to better understand their relationship to UPD proteins from D. melanogaster, we searched for genes encoding UPD homologs in other bee genomes using the three honey bee UPDl amino acid sequences, and constructed a phylogenetic tree. These genomes represent both social and solitary bee species broadly distributed phylogenetically within bees. The LEPTIN protein from Homo sapiens was used as an outgroup (Londraville et al., 2017) (Fig. 4B). As expected, we found that the bee UPDl proteins cluster as three groups that match their location in the UPD locus in the same gene order found for A. mellifera, namely UpdlA, UpdlB, UpdlC (Table S2). Contrary to our expectations, these bee genes did not cluster with the three Drosophila genes, but instead formed a single clade.

Hippo has been implicated in regulating multiple Upd genes in Drosophila (Staley and Irvine, 2010; Shaw et al., 2010; Sun and Irvine, 2011; Ren et al., 2013; Houtz et al., 2017). We examined the sequence around the transcriptional start sites (TSS) of the three Upd-like genes (from −2000 bp relative to the TSS to the predicted start codon) for binding sites for Scalloped, the transcription factor that cooperates with Yorkie in response to Hippo activation (consensus sequence NDGHATNT; Zhang et al., 2008) and found a number of sites (Fig. S4). We did not observe any Activator protein 1 (AP-1) binding sites (TGACTCATA; Chatterjee and Bohmann, 2012) in the proximity of the Upd-like genes. However, we did observe a number of consensus sites for Forkhead box O (FOXO) binding (TKTTYACY; Bai et al., 2013), which can also be activated by the JNK pathway (Biteau et al., 2011), in the Upd-like gene loci (Fig. 3C). Interestingly, we found a canonical heat shock element (HSE, consensus sequence GAANNTTCNNGAA; Birch-Machin et al., 2005; Guertin and Lis, 2010) in the gene loci near the transcriptional start site of the UpdCl gene. Full gene sequences of the promoter regions of the Updl genes with key sites are given in Fig. S4.

A number of putative JAK/STAT target genes showed increased expression in our RNA-Seq dataset (Table S1), including the negative regulators Suppressor of Cytokine Signaling at 44Ac (Socs44A) and Ramshackle, the pathway receptor domeless, the intestinal stem cell (ISC) genes Roundabout 2 (Robo2) (Stine et al., 2014) and Sox21a (Zhai et al., 2017), the immune mediator Thioester-containing protein A (TepA) and the EGF ligand Spitz (Zhou et al., 2013). We confirmed the induction of a number of these by qPCR, including Socs44A, Robo2, Sox21a and Spitz (Fig. 5A). We examined the Sox21a gene to look for putative STAT92E binding sites. In mammals, STAT proteins bind DNA at gamma interferon activation site (GAS) elements that contain the palindromic sequence TTC(n)GAA, where n represents a spacer. The Drosophila Stat92E is able to bind to and activate target gene expression through GAS with either 3 or 4 nucleotide spacers (Rivas et al., 2008). We observed a number of putative sites for STAT92E binding in the Sox21a gene (Fig. 5B).

Alterations in the expression of genes of the JNK and Hippo pathways suggested that heat shock had caused tissue damage that induced a regenerative response. In order to examine this possibility, we first looked at the structure of the midgut after thermal stress. Hematoxylin and Eosin stained sections of midgut tissue revealed a dramatic sloughing of cells into the lumen in the midguts of bees exposed to 45°C for 4 h, which was not observed in bees maintained at 35°C (Fig. 6). Sloughing of damaged cells has been observed in the fruit fly intestine after tissue damage (Buchon et al., 2010). Sloughing had ceased and histological comparison did not reveal any differences between bees maintained at 35 or 45°C for 4 h by 20 h after cessation of thermal stress (data not shown).

To determine whether damage induced by heat shock resulted in cell death in midgut epithelial cells, we used a TUNEL assay to find evidence of DNA breaks occurring as part of the programmed cell death process. We observed increased numbers of cells containing Tdt-labelled nuclei in the midguts of bees maintained at 45°C for 4 h compared with the near absence of these cells in the midguts of bees kept at 35°C for 4 h (Fig. 7). Apoptotic cells were observed in the epithelial layer as well as in the lumen, suggesting that apoptotic (or necrotic) cells are shed concurrent with cell death.

Despite the evidence of significant damage to the midgut caused by exposure to 45°C for 4 h, we did not observe a difference in survival of bees maintained at the two temperatures for up to 8 days after conclusion of the thermal stress (Fig. 1A). While long-term impacts on organismal fitness might be expected after such a stress, these results suggest that the tissue regeneration response is capable of returning the midgut to normal function before detrimental effects on medium-term survival occur.

In D. melanogaster, the midgut maintains homeostasis through the constant turnover of cells, all derived from ISCs (reviewed in Jiang et al., 2016). The proliferation and differentiation program of this cell population to maintain homeostasis in response to damage (Amcheslavsky et al., 2009) or infection (Buchon et al., 2009) is regulated by a number of signals emanating from these cells as well as the surrounding tissues (Nászai et al., 2015; Guo et al., 2016; Rodriguez-Fernandez et al., 2020).

In contrast to the detailed understanding of gut homeostasis elucidated in Drosophila, little is known about the molecular mechanisms regulating normal turnover and healing response in the digestive tracts of other insect Orders, including Hymenoptera. For this Order, which among other families includes the >20,000 species of bees worldwide (Michener, 2000), we know significantly less at the cellular level. Multiple groups indicate that the epithelial layer may be organized into higher-order crypt-like structures in honey bees in contrast to the fly (Jimenez and Gilliam, 1989; 1990; Raes et al., 1994). More recent research has shown that that specific cells at the base of the crypt structures are proliferating cells, also suggesting a spatially organized developmental process in these midguts. In honey bees, homeostatic self-renewal of the digestive tract has been shown. This proliferative activity in the digestive tract is influenced by age, social function, diet, chemical stress and infection (Ward et al., 2008; Willard et al., 2011; Forkpah et al., 2014; Panek et al., 2018). However, we have heretofore lacked a comprehensive cell and molecular map of the midgut under normal and stress situations. As the digestive tract is likely to be the initial point of contact between the honey bee and microbial and chemical agents, a clear understanding of the response of this honey tissue to damage induced by these environmental insults is of paramount importance.

Here, we are able to provide the first view of the pathways involved in gut homeostasis in the honey bee, specifically using the response to thermal stress as a model. We report the dramatic transcriptional remodeling of three important cell stress pathways after damage caused by thermal stress: the Hippo pathway, the JNK pathway and the JAK/STAT pathway, all of which have been implicated in damage responses in the fly digestive tract (Hippo: Shaw et al., 2010; Ren et al., 2010; JNK: Biteau et al., 2008; JAK/STAT: Jiang et al., 2009; Buchon et al., 2009; Lin et al., 2010). We found that 14 genes of the Hippo pathway were upregulated. Among these are a number of signaling components as well as the homologs of multiple target genes characterized in Drosophila, including Kibra, a negative feedback regulator of the Hippo pathway, Myc, a gene that is instrumental in promoting proliferation in cells, and potentially the genes encoding the two Upd-like proteins, cytokines that induce ISCs to induce proliferation. Hippo has a well-characterized role in midgut regeneration in the fly (Staley and Irvine, 2010; Karpowicz et al., 2010; Shaw et al., 2010; Ren et al., 2010; Poernbacher et al., 2012; Li et al., 2014). In this context, it appears to possess both cell autonomous and non-cell autonomous roles. When activated in ISCs, it causes proliferation through upregulation of Myc. When turned on in enterocytes, it causes the secretion of UPD molecules that trigger the JAK/STAT pathway in ISCs, further leading to Myc induction and proliferation (reviewed in Jiang et al., 2016). The JNK pathway is key stress response pathway in many systems that has also been implicated in intestinal regeneration (Biteau et al., 2008). Here, we observed upregulation of eight discrete genes associated with the JNK pathway. Among these are a number of signaling components as well as the homologs of multiple target genes characterized in Drosophila, including Puckered, a negative feedback regulator of the JNK pathway. An interaction between the JNK pathway and heat shock has been described in the fly (Gonda et al., 2012) and worm (Oh et al., 2005). There is evidence that in other contexts the Hippo and JNK pathways engage in significant crosstalk (Meng et al., 2016). For example, both Myc (Ren et al., 2013) and the Upd (Biteau et al., 2008) genes are suggested to be JNK targets as well as Hippo targets in the fly intestine. Examination of the promoters of the two induced Upd-like genes revealed sites that are consistent with potential activation of these genes by either the Hippo or JNK pathways. However, in the absence of genome-wide analysis to define site enrichment, ChIP studies to establish factor binding, and loss-of-function studies to demonstrate whether specific sites are necessary or sufficient for gene induction, these findings represent potential starting points for future studies.

In D. melanogaster, the three members of the Unpaired family (UPD1, UPD2 and UPD3) represent the only ligands that activate DOMELESS, the sole receptor of the JAK/STAT pathway in this species (Harrison et al., 1998; Agaisse et al., 2003; Gilbert et al., 2005; Hombría et al., 2005). Activation of this pathway is thought to be one of the two drivers (along with the EGFR MAPK pathway) of proliferation in ISCs in the adult midgut (Nászai et al., 2015; Guo et al., 2016; Rodriguez-Fernandez et al., 2020). The JAK/STAT pathway is also important in later developmental stages for enteroblasts to transition to enterocytes (Lin et al., 2008; Beebe et al., 2009; Jiang et al., 2009). One of the main targets of the JAK/STAT pathway driving proliferation is the Myc gene (Ren et al., 2013), which we found to be upregulated here. All three UPD proteins have been shown to play both overlapping and distinct roles in regulating midgut homeostasis in Drosophila (Lin et al., 2010; Osman et al., 2012; Zhou et al., 2013). Here, we identified and characterized three UPD-like proteins in bees. The homologs of the three UPD proteins have not previously been identified in bees (Evans et al., 2006; Barribeau et al., 2015; Brutscher et al., 2015). All three putative honey bee UPD homologs share a number of highly conserved motifs with the three Drosophila UPD proteins, which may be important for folding, processing, stability or receptor-binding activity. Despite the existence of these conserved regions, the amino acid sequences are highly divergent, which underlies the historical difficulty in finding UPD orthologs in different insect orders (Loehlin and Werren, 2012). A phylogenetic tree based on the full amino acid sequence showed that the bee UPDl proteins cluster as three groups, strongly supported by bootstrapping analysis that matches their conserved order in the UPD locus. Interestingly, the bee sequences formed a single clade, suggesting that the Upd duplications in the lineage leading to bees occurred after its separation from the last common ancestor with Drosophila. However, future work including comparison of UPDl proteins from a broader range of arthropods may allow for a better understanding of the evolution of UPDl genes. While a significant amount of work has been done examining the structural requirements of the most closely related cytokine of vertebrates, LEPTIN, little is known about the role of specific structural elements of the insect UPD proteins (Londraville et al., 2017). Interestingly, despite very low sequence similarity, the folded structure of LEPTIN and the D. melanogaster UPD2 share a similar structure (Londraville et al., 2017), and LEPTIN can actually trigger the D. melanogaster UPD receptor DOMELESS (Rajan and Perrimon, 2012). Future work will be necessary to functionally and structurally characterize these putative UPDl proteins.

This work leaves a number of other key questions to be addressed in future studies. Most crucial will be the demonstration of an increase in proliferation following the damage induced by thermal stress and experiments showing that a regeneration program is indeed an important component of thermotolerance. In addition, we hope to define the functional role of specific candidate genes, such as the Upd-like genes, in this process by leveraging new genetic tools developed for use in honey bees, such as the recently developed method of symbiont-delivered RNAi (Leonard et al., 2020). Similarly, in the work described here, we have not been able to ascribe the observed transcriptional changes or cellular behaviors, such as cell death, to specific cell types. Currently, we have very little information about the cell populations and developmental hierarchy of the honey bee digestive tract and the future the use of in situ hybridization, in conjunction with immunohistochemistry and histological stains, to map specialized cell populations in the midgut to different niches and regions of the digestive tract tissue with and without stress would be highly valuable.

Understanding the impact of environmental stress caused by climate change (Diffenbaugh et al., 2017) on various organisms is critical (Tomanek, 2010). Exposure to temperatures outside of optimal ranges affects animal function, and ultimately fitness, such that the ability to endure temperatures beyond ideal ranges is likely essential (Stillman, 2019). One framework for understanding organismal responses to stress might use the complimentary ideas of resistance and tolerance borrowed from plant ecology and applied to infectious disease in animals to great effect (Råberg et al., 2007; Ayres and Schneider, 2012; Medzhitov et al., 2012). In that context, resistance is defined as the ability to limit the burden of infectious agent (i.e. the stress itself), while tolerance is defined as the ability to withstand the harm caused by the infectious agent (i.e. the damage caused by the stress). These concepts have recently been more broadly applied to explain mechanisms protecting animals from more diverse stresses (Dillman and Schneider, 2015) and can be further applied to thinking about mechanisms protecting organisms from thermal stress. In the context of thermal stress, thermoresistance would entail any behavioral or physiological mechanisms for maintaining temperature within the optimal range. For individual bees away from the colony, thermoresistance includes the use of liquid regurgitation for evaporative cooling of the head (Heinrich, 1979) and changes in the respiratory rate and abdominal pumping (Kovac et al., 2007). Honey bees (and other social insects) are unusual among insects in that their thermoregulation operates at two levels. While individual honey bees are poikilotherms and their body temperature is greatly influenced by the surrounding air temperature, the colony has homeothermic properties and can maintain a relatively constant temperature under normal conditions. This colony-level homeostatic regulation of hive temperature to maintain it between 32 and 35°C during normal conditions is an important adaptive feature of honey bees (Owens, 1971; Seeley, 1985; Southwick, 1991; Heinrich, 1993; Stabentheiner et al., 2010). Thus, thermoresistance likely takes place largely at the level of the colony, where the limited temperature range is maintained by complex individual behaviors (Stabentheiner et al., 2010, and references therein) and is important for brood development and normal colony function (Himmer, 1927; 1932; Groh et al., 2004; Wang et al., 2016). Climate change has increased the interest in colony-level thermoregulation. For example, a recent study has looked at colony thermoresistance responses to stimulated heat waves and discovered that one adaptive strategy included a significant increase in recruitment of water foragers to increase evaporative cooling (Bordier et al., 2017).

Thermotolerance, in contrast, would include strategies that limit damage caused by thermal stress. At the cellular level, thermotolerance includes HSR-mediated protection of the proteome from thermal stress (Klaips et al., 2018), well characterized in honey bees (McKinstry et al., 2017). However, less is known about mechanisms of thermotolerance operating at the tissue and organismal levels and no information exists on how thermotolerance might operate at the colony level. Here, we describe results that implicate a robust tissue regeneration program in supporting thermotolerance at the tissue and organismal levels. Specifically, thermal stress induces a dramatic induction in signaling pathways associated with regeneration of the digestive tract in the fruit fly. Cell death and shedding of epithelial cells is apparent in the midguts of thermally stressed bees, likely triggering this apparent regenerative program. Despite the evidence of significant damage to the midgut caused by 4 h at 45°C, the midgut had returned to normal by 20 h post-thermal stress and we did not observe a difference in survival of bees maintained at the two temperatures for up to 8 days after conclusion of the thermal stress, a finding which is consistent with published critical thermal maximum values for honey bees (Kovac et al., 2014; Li et al., 2019). While long-term impacts on organismal fitness might be expected after such a stress, these results suggest that the tissue regeneration response is capable of returning the midgut to normal function before detrimental effects on medium-term survival occur. While this study suggests that a strong tissue regeneration program may contribute to thermotolerance at the tissue and organismal level, further work will be important to determine its relative importance in honey bee thermotolerance overall. Furthermore, this work provides significant initial characterization of the damage response and regeneration program in this key pollinating species, with implications for other related species.

The authors thank Stephen Chan for helpful comments and critical review of the manuscript.