\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Tissue-specific transcriptional patterns underlie seasonal phenotypes in
honey bees (Apis mellifera)}



\author[1]{Sean Bresnahan}%
\author[2]{Mehmet Döke}%
\author[2]{Tugrul Giray}%
\author[1]{Christina Grozinger}%
\affil[1]{The Pennsylvania State University}%
\affil[2]{University of Puerto Rico}%


\vspace{-1em}



  
  \date{May 24, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Faced with adverse conditions, such as winter in temperate regions or
hot and dry conditions in tropical regions, many insect species enter a
state of diapause, a period of dormancy associated with a reduction or
arrest of physical activity, development, and reproduction. Changes in
common physiological pathways underlie diapause phenotypes in different
insect species. However, most transcriptomic studies of diapause have
not simultaneously evaluated and compared expression patterns in
different tissues. Honey bees (Apis mellifera) represent a unique model
system to study the mechanisms underpinning diapause. In winter, honey
bees exhibit a classic diapause phenotype, with reduced metabolic
activity, increased physiological nutritional resources, and altered
hormonal profiles. However, winter bees actively heat their colony by
vibrating their wing muscles; thus, this tissue is not quiescent. Here,
we evaluated the transcriptional profiles of flight muscle tissue and
fat body tissue (involved in nutrient storage, metabolism and immune
function) of winter bees. We also evaluated two behavioral phenotypes of
summer bees: nurses, which exhibit high nutritional stores and low
flight activity, and foragers, which exhibit low nutritional stores and
high flight activity. We found winter bees and nurses have similar fat
body transcriptional profiles compared to foragers, whereas winter bees
and foragers have similar flight muscle transcriptional profiles
compared to nurses. Additionally, differentially expressed genes were
enriched in diapause-related GO terms. Thus, honey bees exhibit
tissue-specific transcriptional profiles associated with diapause,
laying the groundwork for future studies evaluating the mechanisms,
evolution, and consequences of this tissue-specific regulation.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Diapause is a genetically programmed state of dormancy, which is
triggered by environmental cues in advance of adverse conditions, such
as winter in temperate areas or hot and dry conditions in tropical areas
(Denlinger, D. L. \& Armbruster, P. A., 2014, Saulich, A. \& Musolin,
D., 2017; Sim, C. \& Denlinger, D. L., 2013). Many insect species
exhibit a diapause phenotype, which typically involves a reduction or
arrest in physical activity, development and/or reproduction. Different
insect species diapause in different developmental stages, from eggs to
adults. Several studies have evaluated the transcriptional and
physiological processes underpinning diapause initiation, maintenance,
and termination in different species (Amsalem, E. et al., 2015; Ragland,
G. J. \& Keep, E., 2017; Ragland, G. et al., 2019; Santos, P. K. F. et
al., 2018; Treanore, E. D. et al., 2020; Yocum, G. \& Rinehart, J.,
2015). Common physiological processes have been found to be associated
with the diapause phenotype, across developmental stages and insect
species. These include changes in regulation of energetic resources
(e.g., increased feeding, reduced metabolism, upregulation of fatty acid
synthesis), altered stress and immune responses, changes in circadian
rhythm, and altered hormonal profiles. Transcriptomic studies of
diapause phenotypes across species have demonstrated that there is no
common set of genes underpinning these processes: rather, species use
different genes, albeit in similar functional categories, to generate a
diapause phenotype (Amsalem, E. et al., 2015; Ko\selectlanguage{polish}š\selectlanguage{english}t\selectlanguage{ngerman}ál, V. et al., 2017;
Ragland, G. J. \& Keep, E., 2017; Treanore, E. D. et al., 2020).
However, the full range of variation in molecular processes underpinning
the diapause phenotype has not been evaluated, since studies largely
focus on whole-body or single tissue expression patterns: within the
same individual, different tissues may span the range from quiescent to
active.

Honey bees (\emph{Apis mellifera} ) are a unique species in which to
study the molecular and physiological pathways underpinning diapause.
Honey bees live in perennial colonies, with a single reproductive female
queen and tens of thousands of facultatively sterile female workers
(Winston, M. L., 1987). Honey bees evolved in tropical regions, and
later expanded into temperate Europe (reviewed in Dogantzis, K. A. \&
Zayed, A., 2019). In both temperate and tropical regions, bees
experience periods of season dearth during which they reduce or cease
brood rearing, requiring that the adult workers exhibit an increased
lifespan. This long-lived, diapause-like phenotype has been best studied
in winter conditions (reviewed in Döke, M. A. et al. 2015 and Grozinger,
C.M. et al., 2014). During the temperate winter, honey bees feed on
stored honey and pollen and form a thermoregulating cluster. Colonies
cease brood rearing in the winter, and winter worker bees exhibit a
unique physiological phenotype and significantly extended lifespan
compared to summer bees. In tropical regions, honey bees will cease
brood production during hot and dry conditions, and can leave their
hives to migrate long distances before establishing a new colony
(Grozinger, C. M., et al., 2014). Change in worker longevity has not
been studied in these migrating colonies. However, in the absence of
young bees to replace the older bees, the existing adult workers
presumably must live longer for the colony to survive. Indeed, our
recent studies of honey bees on the tropical island of Puerto Rico
demonstrated that worker bees' lifespan changes throughout the year and
long-lived worker bee phenotype can be induced by removing the brood
from colonies at any time in the year (Feliciano-Cardona, S. et al.,
2020). Other studies have demonstrated that caging the queen and
preventing her from laying eggs can result in long-lived worker bees,
even in summer months in temperate climates (reviewed in Döke, M. A. et
al., 2015).

The physiology of winter worker bees (and long-lived summer bees)
suggest that they are in a diapause state. Winter worker bees will live
up to eight months, whereas summer worker bees typically live for only a
few (\textasciitilde{}6) weeks (Fluri et al. 1982). In the summer
months, worker bees will transition between different physiological and
behavioral states (reviewed in Robinson, G. E., 1992). When they are
young, honey bees act as nurses, feeding developing brood from
secretions produced by specialized glands in their heads. When they are
middle-aged (\textasciitilde{}2 weeks old), worker bees will transition
to other tasks in the colony, such as comb-building or guarding. In
their final behavioral phase, worker bees will transition to become
foragers, where they collect resources from the surrounding landscape.
Nurse bees are characterized by having large nutritional (fat) stores in
their abdominal fat bodies, and the loss of these stores is associated
with (and can accelerate) the transition to foraging (Toth, A. L. \&
Robinson, G. E., 2005; Toth, A. L. et al., 2005). Compared to foragers,
nurses also have significantly higher hemolymph levels of vitellogenin
(Vg), a nutritional storage protein that is converted into brood food
(Amdam, G. V. et al., 2003). Vg also negatively regulates levels of a
key hormone, juvenile hormone (JH), and JH levels are significantly
lower in nurse bees than in forager bees (reviewed in Amdam, G. V. et
al., 2009). Decreasing levels of Vg or increasing levels of JH can
accelerate the transition to foraging behavior in worker bees. Like
summer nurse bees, long-lived winter honey bees exhibit high nutritional
stores, high levels of Vg, and low levels of JH (reviewed in Döke, M. A.
et al., 2015).

However, winter bees also exhibit similarities to summer foragers, in
terms of their use of their thoracic flight muscles. Winter bees
actively use their wing muscles to generate heat when exterior
temperatures drop below 10°C; this process is so effective that worker
bees can increase the internal temperatures of the hive to 34°C in late
winter/early spring to support brood rearing (Döke, M. A. et al., 2015;
Seeley, T. D. \& Visscher, P. K., 1985; Currie, R. W. et al., 2015).
Forager bees are active fliers and can travel several kilometers from
the colony during foraging trips (Couvillon, M. J. et al., 2015). Thus,
it could be hypothesized that winter bees and foragers may have flight
muscles that have distinct transcriptional, physiological or biochemical
profiles than nurse bees. Indeed, studies have found differences in the
transcriptional profiles of nurse and forager bees is primarily affected
by age, and not by behavioral state or even foraging experience
(Margotta, J. W. et al., 2012; Schippers, M-P. et al., 2010). The
improved flight ability of foragers may be because they have lost
considerable mass during the transition to foraging, by reducing their
nutritional stores (Vance, J. T. et al., 2009). However, there are
differences in oxidative damage which are generated by flight muscle use
(Margotta, J. W. et al., 2018) which may be associated with differences
in expression patterns in flight muscles based on activity (Oskay, D.,
2007).

The genome-wide transcriptional profile associated with the winter bee
physiological state has not been characterized. One challenge with such
studies is that it is unclear which ``physiological state'' can serve as
a non-diapausing comparison, since summer worker bees have two distinct
physiological states which vary significantly in physiological processes
classically associated with diapause (nurses versus foragers). It could
be expected that winter bees and nurse bees would resemble each other,
as states with considerable nutritional resources and the ability to
enter a long-lived phenotype. However, winter bees and forager bees are
both capable of extended use of their flight muscles, while nurses
typically are not active fliers. Interestingly, \emph{k} -means
clustering analysis of candidate gene expression profiles (n = 9 in fat
body tissues and n = 5 in flight muscle tissues) demonstrated that
winter bee fat body fat body tissue profiles clustered with nurses,
while winter bee flight muscle tissue profiles clustered with foragers
(Döke, M. A. et al., in preparation). Similarly, brain expression
patterns of long-lived summer bees collected from brood-less colonies
resembled nurse bees more than foragers (Münch, D. et al., 2013;
Whitfield, C. et al., 2003). Overall, these studies would suggest that
there are expression profiles in winter bees that are consistent with
physiological processes that differ between summer nurses and foragers,
some of which are classically associated with the diapause phenotype.
However, these results also suggest there is considerable variation
across tissues.

Here, we evaluated and compared the transcriptional profiles of winter
bees (collected in December), summer nurse and forager bees (collected
in July), in abdominal samples containing fat body tissue and thoracic
samples containing flight muscles. First, we determined the extent to
which these two tissues varied in expression across these different
phenotypes. Based on previous studies, we expected large differences in
fat body expression patterns and much fewer differences in flight muscle
tissue. Second, we determined if and how the expression patterns of
winter bees corresponded to those of summer bees, predicting that winter
bee fat body profiles will be more similar to summer nurse fat body
tissue profiles versus foragers, while winter bee flight muscle tissue
will be more similar to foragers. Third, we determined if the fat body
and flight muscle transcriptional profiles showed regulation of
functional categories of genes previously associated with diapause in
different insects species -- e.g., metabolic processes, stress response,
circadian rhythm (Amsalem, E. et al., 2015; Denlinger, D. L., \&
Armbruster, P. A., 2014; Ko\selectlanguage{polish}š\selectlanguage{english}t\selectlanguage{ngerman}ál, \selectlanguage{polish}Š, \selectlanguage{english}T. et al., 2017; Kunk, M. et al.,
2019; Ragland, G. J. \& Keep, E., 2017; Santos, P, K. F. et al., 2018;
Treanore, E. D. et al., 2020; Yamada, H. \& Yamamoto, M. T., 2011;
Yocum, G. D. et al., 2018). We predicted the fat body tissues will show
classic signatures of diapause, but the flight muscle tissue will not.

\textbf{Methods}

\emph{Sample collection and RNA-sequencing}

Colonies for this study were generated and described as in D\selectlanguage{ngerman}öke M. A. et
al. (2019). Sample collection was conducted as described in Döke, M. A.
(2017). Briefly, four colonies headed by naturally mated honey bee
queens from four different commercial breeders were established at a
Pennsylvania State University apiary in central Pennsylvania. Colonies
were maintained according to standard apicultural practices. More than 6
weeks after establishment (July 2013), which is sufficient time for the
eggs laid by the queens to develop into adult nurses and foragers, 6
summer nurse and forager bees per colony were collected on dry ice and
stored at -80\textsuperscript{o}C. Bees were collected based on
behavioral observations, with nurses collected as they were observed
feeding young brood, while foragers were collected at the colony
entrance as they returned with pollen loads. Winter bees were collected
5 months later (December 2013) from within their winter clusters.

After removal from storage at -80\textsuperscript{o}C, sampled bees were
submerged in Invitrogen\textsuperscript{TM}
Ambion\textsuperscript{TM}RNA\emph{later} \textsuperscript{TM}
Stabilization Solution (ThermoFisher Scientific) to prevent RNA
degradation, and dissected on a platform surrounded by dry ice. To
collect flight muscle tissue, legs and wings were detached from the
thorax and discarded, leaving flight muscles and exoskeleton. To collect
fat body tissues, abdomens were allowed to thaw RNA\emph{later}
\textsuperscript{TM} on ice. Digestive tracts and reproductive organs
were removed by gently pulling from the stinger using forceps, leaving
the eviscerated abdomen with fat bodies attached to the exoskeleton.
Dissected tissues for each worker were individually placed in a 1.5ml
nuclease-free microcentrifuge tube (Denville Scientific Inc., Metuchen,
NJ) and kept at -80\textsuperscript{o}C.

For each of the 4 colonies, samples of the same tissue and phenotype
were pooled (6 bees per sample). Pooled samples (4 colonies x 3 groups x
2 tissues = 24 samples, total) were homogenized using an automated
homogenizer (Thermo Savant FastPrep\textsuperscript{TM} FP120 Cell
Disrupter System) for 3 intervals of 45 seconds at the highest speed
setting, using 3-6 2mm zirconia beads (BioSpec Products, Inc.).
Homogenates were stored for later use in RNA extraction via the RNeasy
Plus Universal kit (QIAGEN) for RNA sequencing (RNA-seq). RNA libraries
were prepared by the Sequencing and Genomic Technologies Shared Resource
at the Duke Center for Genomic and Computational Biology, Durham, NC,
and sequenced on an Illumina HiSeq 4000 platform (see Table S1 for a
list of samples).

\emph{Differential gene expression analysis}

NGS reads were assessed for quality control metrics after adapter
trimming with Trim Galore (Krueger, F., 2019). We then used Kallisto
(Bray, N. L. et al., 2016) to quantify transcript abundance by
pseudo-alignment on the latest release of the \emph{Apis
mellifera}genome assembly Amel\_HAv3.1 (Wallberg, A. et al., 2019),
compiled the transcript abundance estimates for each sample into an
expression matrix in R (R Core Team, 2019), and summed transcripts by
gene using the tximport method, see Table S2 (Soneson, C. et al., 2015).
Low count genes (\textless{} 1 count across \textgreater{} 25\% of
samples) were removed from the analysis, resulting in 92.38\%
(9178/9935) of the\emph{A. mellifera} protein-coding genes (Wallberg, A.
et al., 2019) represented in this study.

We employed a two-factor design using DESeq2 (Love, M. I. et al., 2014)
to assess differential expression of protein-coding genes (DEGs) within
tissue types between nurse, forager, and winter bees, using colony as a
cofactor. In each analysis, we applied the default Benjamini-Hochberg
(Benjamini, Y. \& Hochberg, Y., 1995) correction and set the threshold
for differential expression at \emph{padj} \textless{} 0.01.

To assess sample relationships and variance, hierarchical clustering
analysis was conducted using the variance stabilized transformed (VST)
counts generated by DESeq2. We determined that two of the nurse fat body
samples (N1 and N3) were outliers from this analysis, as each created
separate branches independent of the other samples that clustered as
expected by phenotype (Figure S1). We identified 411 genes that were
differentially expressed between the clustered and outlier nurse samples
(\emph{padj} \textless{} 0.05) but did not observe functional enrichment
for any biological processes among these genes. FastQC did not report
major differences in any of the default sequence quality metrics apart
from reads corresponding to highly overrepresented sequences. A local
BLASTN (McGinnis, S. \& Madden, T. L., 2004) search of the
overrepresented sequences revealed many hits with 100\% identity and
100\% query cover for \emph{A. mellifera} 16S ribosomal RNAs, a common
source of contamination during library preparation (Zhao, W. et al.,
2014). However, upon filtering our sequencing reads for honey bee rRNAs
(The RNAcentral Consortium, 2018), sample clustering did not improve
(Figure S2). Thus, we could not determine the true source of this
discrepancy, and so removed these samples from all downstream analyses.

\(\chi^{2}\) \emph{tests and k-means clustering}

We performed chi-squared tests (with Yate's continuity correction) to
determine, for each tissue, whether the number of DEGs between winter
bees and nurses is significantly different than the number of DEGs
between winter bees and foragers. Additionally, we performed\emph{k}
-means clustering (Oyelade, J. et al., 2016) to assess whether the
tissue-specific transcriptomes of winter bees had more in common with
foragers or nurses using the fviz\_cluster function from the factoextra
(Kassambara, A., 2019) package in R, testing for
\(k=2\)clusters. Separate tests were performed using: 1) all
genes passing our initial low-count filter (9178), and 2) only forager
vs nurse DEGs.

\emph{Gene ontology enrichment analyses}

Gene ontology enrichment analyses were performed with topGO (Alexa, A.
et al., 2006), using \emph{Drosophila melanogaster} orthologs of honey
bee genes converted via a one-to-one reciprocal best hit BLAST (Table
S15). We focused our assessment (Tables 3-4) on the more specific and
highly enriched terms through sorting each resulting GO term list by
percent coverage (defined as the proportion of genes associated with a
GO term represented among the differentially expressed genes in each
phenotype contrast) and selecting the top quartile, and then by GO term
rank (defined as the number of child terms per GO term) and selecting
the bottom quartile. Thus, Tables 3-4 represent the significantly
enriched GO categories that had the largest proportion of DEGs and the
most specific biological functions.

\emph{Conservation of biological processes associated with diapause
across insect species}

GO terms found to be associated with diapause in bumble bees
(\emph{Bombus terrestris} ) (Amsalem, E. et al., 2015), an
oil-collecting bee \emph{Tetrapedia diversipes} (Santos, P. K. F. et
al., 2018), alfalfa leaf-cutting bees (\emph{Megachile rotundata} )
(Yocum, G. D. et al., 2018), and across 11 insect species spanning
Diptera, Lepidoptera, and Hymenoptera (Ragland, G. J. \& Keep, E., 2017)
were compiled from previous studies (Table S16). We then tested whether
the enriched GO terms in our study overlapped with these terms.

\textbf{Results}

\emph{Evaluation of the summer and winter bee transcriptomes}

Through RNA-sequencing we obtained on average 17.3 million
pseudo-alignments on the honey bee reference transcriptome for each
biological replicate of nurse (N), forager (F) and winter (W). For the
protein coding genes in the honey bee genome (Amel\_HAv3.1) (Wallberg,
A. et al., 2019), 92.38\% (9178/9935) passed our low count filter.

The number of genes significantly differentially expressed (DEGs)
between nurse, forager, and winter bees in each tissue are listed in
Table 1. In fat body tissue, we found 683 genes differentially expressed
between nurses and foragers, 1617 between winter bees and foragers, and
934 between winter bees and nurses. In flight muscle flight muscle
tissue, we found 2673 DEGs between nurses and foragers, 1176 between
winter bees and foragers, and 2855 between winter bees and nurses.
Interestingly, there were many more DEGs in flight muscle tissues (2763)
versus abdominal tissue (693) in nurses versus foragers. (Table S3-S7)

We performed chi-squared tests to determine, for each tissue, whether
the number of DEGs between winter bees and nurses is significantly
different than the number of DEGs between winter bees and foragers. The
results are displayed in Table 2. There are significantly more genes
differentially expressed between winter bees and forager bees (versus
nurse bees) in the fat body tissues
(\(\chi^{2}=211.76\),\(p<{2.26*10}^{-16}\)), suggesting that winter bees are
more similar to nurse bees in this tissue. In the flight muscle, there
are significantly more genes differentially expressed between winter
bees and nurse bees (relative to forager bees;
\(\chi^{2}=895.06\),\(p<{2.26*10}^{-16}\)), suggesting that winter bees are
more similar to forager bees in this tissue.

To further evaluate how the expression patterns of winter bees
correspond to those of summer bees, we conducted hierarchical clustering
(Figure 1) and \emph{k} -means clustering analyses (calling k = 2
clusters), using all 9178 genes (Figure 2). Both analyses demonstrated
that, in fat body tissue, winter bee and nurse samples form a cluster
independently of foragers, while in flight muscle tissue, winter bee and
forager samples form a cluster independently of nurses. When clustering
with only genes that were differentially expressed between nurses and
foragers (n = 2763 flight muscle tissue DEGs; n = 683 fat body tissue
DEGs) (Figure S3), the same clustering is observed.

\emph{Evaluation of functional categories of summer and winter bee
differentially expressed genes}

Gene ontology (GO) enrichment analyses showed that several biological
processes are differentially regulated between each group (Tables
S9-14). The most significantly enriched and specific terms for each
comparison are displayed in Tables 3-4.

In fat body tissues, DEGs between nurses and foragers were enriched for
GO terms associated with the regulation of dendrite development and
morphogenesis, mRNA processing, and feeding behavior. DEGs between
winter bees and foragers were enriched for GO terms associated with
adult feeding behavior, thermotaxis, and neuromuscular junction
development, in addition to transcriptional regulatory processes,
including mRNA catabolism and poly(A) tail shortening. In contrast, DEGs
between winter bees and nurses were enriched for terms associated with
lumen formation of open tracheal systems, actin filament capping, muscle
filament assembly, and germarium-derived oocyte differentiation and fate
determination (Table 3).

In flight muscles, DEGs between nurses and foragers were enriched for GO
terms associated with regulation of humoral immune response, somatic
muscle development, RNA stability, and membrane protein organization.
DEGs between winter bees and foragers were enriched for GO terms
associated with protein maturation by iron-sulfur cluster transfer,
regulation of rhodopsin mediated signaling pathway, and regulation of
compound eye retinal cell death. DEGs between winter bees and nurses
were enriched for GO terms associated with several processes related to
morphological development, protein localization, post-translational
protein modifications, and immune system function (Table 4).

\emph{Conservation of biological processes associated with diapause
across insect species}

GO terms found to be associated with diapause in bumble bees
(\emph{Bombus terrestris} ) (Amsalem, E. et al., 2015), an
oil-collecting bee \emph{Tetrapedia diversipes} (Santos, P. K. F. et
al., 2018), alfalfa leaf-cutting bees (\emph{Megachile rotundata} )
(Yocum, G. D. et al., 2018), and across 11 insect species spanning
Diptera, Lepidoptera, and Hymenoptera (Ragland, G. J. \& Keep, E., 2017)
were compiled from previous studies. We compared these terms with those
associated with DEGs between summer and winter bees in our study. We
found that one category from fat bodies, dendrite morphogenesis
(GO:0048813), overlapped, and two categories, neurogenesis (GO:0022008)
and morphogenesis of an epithelium (GO:0002009), from flight muscle
tissues overlapped with previous studies.

Among genes related to dendrite morphogenesis, 31 were differentially
regulated between winter and summer bees in the fat bodies, while 39
were differentially regulated in the flight muscles (though note that
this was not a significantly enriched GO category in this tissue). In
total, 49 genes associated with dendrite morphogenesis were
differentially regulated across both tissues between winter and summer
bees, including homologs of genes encoding TGF-beta signaling pathway
components DAF-1, DAF-4, and DAF-8 (TGF-beta receptor type 1
{[}\emph{LOC550930} {]}, activin receptor type 2A {[}\emph{LOC412471}
{]}, and mothers against decapentaplegic homolog 3 {[}\emph{LOC412601}
{]}).

Among genes related to neurogenesis, 3 were differentially regulated
between winter and summer bees in the fat bodies, while 4 were
differentially regulated in the flight muscles. In total, 6 genes
associated with neurogenesis were differentially regulated across both
tissues between winter and summer bees, including epidermal growth
factor receptor-like {[}\emph{LOC100577393} {]}, histone-lysine
N-methyltransferase E(z) {[}\emph{LOC552235} {]}, sn1-specific
diacylglycerol lipase alpha {[}\emph{LOC413224} {]}, protein
gooseberry-neuro-like {[}\emph{LOC411374} {]}, BTB/POZ domain-containing
protein 6-B {[}\emph{LOC100578209} {]}, and putative
glutathione-specific gamma-glutamylcyclotransferase 2
{[}\emph{LOC408369} {]}.

Lastly, 6 genes related to morphogenesis of an epithelium were
differentially regulated between winter and summer bees in the fat
bodies, while 5 were differentially regulated in the flight muscles. In
total, 6 genes associated with morphogenesis of an epithelium were
differentially regulated across both tissues between winter and summer
bees, including epidermal growth factor receptor-like
{[}\emph{LOC100577393} {]}, cyclin-dependent kinase 14
{[}\emph{LOC552617} {]}, radixin homolog 1 {[}\emph{LOC412799} {]},
mitochondrial electron transfer flavoprotein subunit alpha
{[}\emph{LOC551710} {]}, and dystroglycan {[}\emph{LOC408826} {]}.

\textbf{Discussion}

Our results demonstrate that winter bees and nurses have similar
transcriptional profiles in fat body tissues relative to foragers, while
winter bees and foragers have similar transcriptional profiles relative
to nurses. Thus, there are clearly tissue-specific expression patterns
associated with the winter bee phenotype. These patterns are consistent
with the differential functions of these tissues in bees. Winter bees
thus represent a ``mix and match'' phenotype between summer nurse bees
and summer forager bees, with winter fat bodies serving to store
nutrition (as is the case of nurse bees) and winter flight muscles
remaining active (as in forager bees) to generate heat. Our analysis of
biological processes demonstrated that genes differentially expressed
within both tissues corresponded to categories associated with the
diapause phenotype in other insect species.

We found that there were twice as many differentially expressed genes in
flight muscle tissue than in fat body tissues (Table 1). While many
studies have evaluated differences between nurses and foragers in fat
body tissue (Ament, S. A. et al., 2012; Khamis, A. M. et al., 2015;
Seehuus, S-C. et al., 2013), few have evaluated flight muscle tissue.
Interestingly, studies of flight muscle tissue report large changes in
muscle performance or physiology occurring in young, nurse-age, worker
bees, which are more related to maturation processes than the transition
from nursing to forager (Harrison, J., 1986; Herold, R., Borei, H.,
1963; Schippers, M-P. et al., 2010). This would suggest that few genes
would vary between older nurses and foragers, which were used in our
study. However, we found hundreds of differentially expressed genes in
the flight muscle of nurses, foragers and winter bees. In our samples,
bees were actively performing the behaviors (foragers were collected as
they returned to the hive, and winter bees were collected from
thermoregulating clusters) and thus expression patterns may be showing
differences between active versus quiescent muscle tissue. However,
overall transcriptional profiles of muscle tissue have not been broadly
investigated, and future studies should evaluate to what extent there
are ``baseline'' differences in expression patterns in flight muscle
tissue in these different phenotypes, which may prime the bees for
different levels of activity.

Based on total numbers of DEGs and clustering analyses, our results
suggest that winter bees' transcriptional profiles in the fat bodies
that support the long-lived winter physiological state, while
transcriptional profiles in the flight tissues support (or are
reflective of) high levels of activity. Thus, winter bees appear to be
``mixing and matching'' the gene expression profiles and underlying
physiological processes of nurse and forager bees, in a tissue-specific
way. These results highlight the importance of considering
tissue-specific expression patterns when evaluating processes involved
in diapause. It will be valuable to determine how this tissue-specific
phenotypic plasticity is controlled, both at the molecular and
behavioral level. Nurse-like transcription profiles in the fat body
tissue may be a result of reduced exposure to brood pheromone (since the
there are no developing larvae in the colony) and potentially the
presence of older forager-age bees which release ethyl oleate. Brood
pheromone can accelerate the transition from nursing to foraging, reduce
Vitellogenin levels and reduce longevity in worker bees (Amdam, G. V. et
al., 2009). Exposure to ethyl oleate can slow the behavioral maturation
from nursing to forager (Leoncini, I. et al., 2004). Forager-like
transcriptional profiles in the muscle tissue may be the result of
muscle activity, but it remains to be determined what neurophysiological
mechanism triggers the generation of muscle activity in these winter
bees.

When comparing the GO categories identified in our study with those from
previous studies, we found differential regulation of genes associated
with dendrite morphogenesis (GO:0048813) neurogenesis (GO:0022008) and
morphogenesis of an epithelium (GO:0002009). Several of the
differentially regulated dendrite morphogenesis-related genes
participate in mTOR and TGF-beta signaling. Interestingly, dendrite
morphological restructuring, driven in part by TGF-beta signaling,
underlies a related occurrence of phenotypic plasticity in \emph{C.
elegans} : the Dauer phenotype. The Dauer larva is a stress-induced
developmental stage wherein several tissues (particularly nervous
tissue) undergo remodeling, producing a stress-resistant and long-lived
alternative developmental phenotype, akin to the diapause state
(Androwski, A. et al., 2017). Molecular analyses have established the
Dauer stage to be driven by genes encoding components of an
insulin-related pathway, a cyclic nucleotide pathway, and a
TGF-beta-related pathway. DAF-7 -- Dauer larva development regulatory
growth factor 7 -- is the TGF-beta-related ligand for the Dauer pathway.
The DAF-7 signal is transduced by the DAF-1 and DAF-4 TGF-beta-family
type 1 receptors as well as several SMAD-family transcription factors,
including DAF-8. The honey bee orthologs of genes encoding DAF-1, DAF-4,
and DAF-8 are TGF-beta receptor type 1 {[}\emph{LOC550930} {]}, activin
receptor type 2A {[}\emph{LOC412471} {]}, and mothers against
decapentaplegic homolog 3 {[}\emph{LOC412601} {]}, which were
differentially regulated between winter bees and foragers in fat bodies.
Similarly, the Dauer phenotype has been shown to promote neurogenesis of
mechanosensory neurons that were experimentally damaged in \emph{C.
elegans} (Caneo, M. et al., 2019). Interestingly, studies in the
vertebrate African turquoise killifish (\emph{Nothobranchius furzeri} )
(Chi-Kuo H. et al., 2020) and invertebrate cotton bollworm
(\emph{Helicoverpa armigera} ) (Yu-Xuan L. et al., 2013), species which
both exhibit diapause phenotypes, reveal a role for
neurogenesis-associated Polycomb repressive complex members such as
histone-lysine N-methyltransferase E(z) {[}\emph{LOC552235} {]} (which
was upregulated in winter bees and in foragers compared to nurses in the
flight muscles) in diapause maintenance. Thus, conserved signaling and
regulatory pathways might underly stress-induced long-lived
developmental phenotypes across taxa.

Overall, our results indicate that honey bees exhibit tissue-specific
transcriptional profiles associated with diapause as an adaptation to
different seasonal conditions. Differential regulation of genes
associated with dendrite morphogenesis and neurogenesis, including
members of the TGF-beta signaling pathway and the Polycomb repressive
complex suggest that conserved molecular pathways may underlie
stress-induced long-lived developmental phenotypes across taxa. These
studies lay the groundwork for future evaluations of the mechanisms,
evolution, and consequences of these interrelated phenomena.

\textbf{Acknowledgments}

We would like to thank Maryann Frazier for helping generate the original
colonies used in these studies, and Ryan Reynolds, Sara Ashcraft,
Bernardo Niño, and Mario Padilla for expert beekeeping assistance, and
members of both the Grozinger and Giray labs for feedback on our
manuscript.

\textbf{Conflict of Interest}

The authors declare no conflict of interest associated with any aspect
of this study.

\textbf{References}

Alexa, A., Rahnenführer, J., Lengauer, T. (2006). Improved scoring of
functional groups from gene expression data by decorrelating GO graph
structure. \emph{Bioinformatics.} 22, 1600--1607

Ament, S. A., Wang, Y., Chen C-C., Blatti, C. A., Hong, F. et al.
(2012). The transcription factor \emph{ultraspiracle} influences honey
bee social behavior and behavior-related gene expression. \emph{PLoS
Genetics} . 8(3): e1002596

Amdam, G. V., Norberg, K., Hagen, A., Omholt, S. W. (2003). Social
exploitation of vitellogenin. \emph{PNAS} . 100(4): 1799-1802

Amdam, G. V., Ihle, K. E., Page, R. (2009). Regulation of honeybee
worker (\emph{Apis mellifera} ) life histories by
vitellogenin.\emph{Hormones, Brain and Behavior.} Second Edition.
1003-1027

Amdam, G. V., Rueppel, O., Fondrk, M. K., Page, R. E., Nelson, C. M.
(2009). The nurse's load. Early-life exposure to brood rearing affects
behavior and lifespan in honey bees. \emph{Experimental Gerontology} .
44(6-7): 467-471

Amsalem, E., Galbraith, D. A., Cnaani, J., Teal, P. E. A., Grozinger, C.
M. (2015). Conservation and modification of genetic and physiological
toolkits underpinning diapause in bumble bee queens. \emph{Molecular
Ecology,} 24(22): 5596-5615

Andrews, S. (2019). \emph{FastQC.}

Androwski, R. J., Flatt, K. M., Schroeder, N. E. (2017). Phenotypic
plasticity and remodeling in the stress-induced \emph{C. elegans}
dauer.\emph{Wiley Interdisciplinary Reviews of Developmental Biology} .
6(5): doi:10.1002/wdev.278

Benjamini, Y. \& Hochberg, Y. (1995). Controlling the False Discovery
Rate: A Practical and Powerful Approach to Multiple Testing. \emph{J. R.
Stat. Soc. Ser. B Methodol.} 57, 289--300

Caneo, M., Julian, V., Byrne, A. B., Alkema, M. J., Calixto, A. (2019).
Diapause induces functional axonal regeneration after necrotic insult
in\emph{C. elegans} . \emph{PLoS Genetics} . 15(1): e1007863

Chi-Kuo Hu et al. (2020). Vertebrate diapause preserves organisms long
term through Polycomb complex members. \emph{Science} . 367(6480):
870-874

Currie, R. W., Spivak, M., Reuter, G. S. (2015). Wintering Management of
Honey Bee Colonies. in The Hive and the Honey Bee. \emph{Dadant \& Sons}
.

Couvillon, M.J., Riddell Pearce, F.C., Accleton, C. \emph{et al} .
(2015). Honey bee foraging distance depends on month and forage
type.\emph{Apidologie.} 46, 61--70

Dogantzis, K. A. \& Zayed, A. Recent advances in population and
quantitative genomics of honey bees. \emph{Curr. Op. Insect Sci} . 31,
93--98 (2019)

Denlinger, D. L. \& Armbruster, P. A. (2014). Mosquito
diapause.\emph{Annu. Rev. Entomol.} 59, 73--93

Döke, M. A., Frazier, M., Grozinger, C. M. (2015). Overwintering Honey
Bees: Biology and Management. \emph{Curr. Op. Insect Sci.} 185--193

Döke, M. A. (2017). Ecophysiology of Overwintering in Honey Bees. PhD
thesis.

Döke, M. A., McGrady, C. M., Otieno, M., Grozinger, C. M., Frazier, M.
(2019). Colony size, rather than geographic origin of stocks, predicts
overwintering success in honey bees (Hymenoptera: Apidae) in the
Northeastern United States. \emph{Journal of Economic Entomology} .
112(2): 525-533

Elsik, C. et al\emph{.} (2016). Hymenoptera Genome Database: integrating
genome annotations in HymenopteraMine. \emph{Nucleic Acids Res.} 4,
D793-800

Feliciano-Cardona, Stephanie, \emph{et al} . (2020). Honey bees in the
tropics show winter bee-like longevity in response to seasonal dearth
and brood reduction. \emph{Front. Eco. \& Evo.} 8: 336

Fluri, P., Lüscher, M., Wille, H., Gerig, L. (1982). Changes in weight
of the pharyngeal gland and haemolymph titres of jjuvenile hormone,
protein and vitellogenin in worker honey bees. \emph{Journal of Insect
Physiology} . 28(1): 61-68

Grozinger, C.M., Richards, J., Mattila, H.R. (2014). From molecules to
societies: mechanisms regulating swarming behavior in honey bees (Apis
spp.). \emph{Apidologie.} 45, 327--346

Grozinger, C.M., Richards, J. \& Mattila, H.R. From molecules to
societies: mechanisms regulating swarming behavior in honey bees (Apis
spp.). Apidologie 45, 327--346 (2014).

Harrison, J. (1986). Caste specific changes in honeybee flight
capacity.\emph{Physiological Zoology} . 59: 175-187

Herold, R., Borei, H. (1963). Cytochrome changes during honey bee flight
muscle development. \emph{Developmental Biology} . 8: 67-79

Kassambara, A. (2019). \emph{Factoextra: extract and visualize the
results of multivariate data analyses} .

Khamis, A. M., Hamilton, A. R., Medvedeva, Y. A., Alam, T., Alam, I., et
al. (2015). Insights into the Transcriptional Architecture of Behavioral
Plasticity in the Honey Bee \emph{Apis mellifera} . \emph{Scientific
Reports} , 5: 11136

Ko\selectlanguage{polish}š\selectlanguage{english}t\selectlanguage{ngerman}ál, V., \selectlanguage{polish}Š\selectlanguage{english}t\selectlanguage{polish}ě\selectlanguage{english}tina, T., Poupardin, R., Korbelov\selectlanguage{ngerman}á, J., Bruce, A. W.
(2017). Conceptual framework of the eco-physiological phases of insect
diapause development justified by transcriptomic profiling. \emph{Proc.
Natl. Acad. Sci.} 114, 8532--8537

Krueger, F. (2019). \emph{Trim Galore} .

Kunc, M. et al\emph{.} (2019). The Year of the Honey Bee (Apis mellifera
L.) with Respect to Its Physiology and Immunity: A Search for
Biochemical Markers of Longevity. \emph{Insects.} 10, 244

Leoncini, I., Le Conte, Y., Costagliola, G., Plettner, E., Toth, A. et
al. (2004). Regulation of behavioral maturation by a primer pheromone
produced by adult worker honey bees. \emph{PNAS} . 101(50): 17559-17564

Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold
change and dispersion for RNA-seq data with DESeq2. \emph{Genome
Biol.}15, 550

Margotta, J. W., Mancinelli, G. E., Benito, A. A., Ammons, A., Roberts,
S. P., Elekonich, M. M. (2012). Effects of flight on gene expression and
aging in the honey bee brain and flight muscle. \emph{Insects} . 4(1):
9-30

Margotta, J. W., Roberts, S. P., Elekonich, M. M. (2018). Effects of
flight activity and age on oxidative damage in the honey bee, \emph{Apis
mellifera} . \emph{J. Exp. Bio} . 221, jeb183228

Martin, M. (2011). Cutadapt removes adapter sequences from
high-throughput sequencing reads. \emph{EMBnet.journal.} 17, 10--12

McGinnis, S. \& Madden, T. L. (2004). BLAST: at the core of a powerful
and diverse set of sequence analysis tools. \emph{Nucleic Acids Res.}32,
W20--W25

Münch, D., Kreibich, C. D., Amdam, G. V. (2013). Aging and its
modulation in a long-lived worker caste of the honey bee. \emph{J. Exp.
Bio} . 216(9): 1638-1649

Oskay, D. (2007). Plasticity in flight muscle development and honey bee
division of labor. PhD thesis.

Oyelade, J. et al\emph{.} (2016). Clustering algorithms: their
application to gene expression data. \emph{Bioinform. Biol.
Insights.}10, 237--253

R Core Team (2019). \emph{R: A language and environment for statistical
computing} . (R Foundation for Statistical Computing)

Ragland, G. J., Armbruster, P. A., Meuti, M. E. (2019). Evolutionary and
functional genetics of insect diapause: a call for greater
integration.\emph{Curr. Op. Insect Sci.} (2019). 36: 74-81

Ragland, G. J. \& Keep, E. (2017). Comparative transcriptomics support
evolutionary convergence of diapause responses across Insecta:
Comparative transcriptomics of diapause. \emph{Physiol. Entomol.} 42,
246--256

Robinson, G. E. (1992). Regulation of division of labor in insect
societies. \emph{Annu. Rev. Entomol.} 37, 637--665

Santos, P. K. F., de Souza Araujo, N., Françoso, E., Zuntini, A. R.,
Arias, M. C. (2018). Diapause in a tropical oil-collecting bee:
molecular basis unveiled by RNA-Seq. \emph{BMC Genomics} 19, 305

Saulich, A. Kh. \& Musolin, D. L. (2017). Summer diapause as a special
seasonal adaptation in insects: diversity of forms, control mechanisms,
and ecological importance. \emph{Entomological Review} . 97(9):
1183-1212

Seehuus, S-C., Taylor, S., Petersen, K., Aamodt, R. M. (2013). Somatic
maintenance resources in the honeybee worker fat body are distributed to
withstand the most life-threatening challenges at each life
stage.\emph{PLoS One} . 8(8): e69870

Seeley, T. D. \& Visscher, P. K. (1985). Survival of honeybees in cold
climates: the critical timing of colony growth and
reproduction.\emph{Ecol. Entomol.} 10, 81--88

Schippers, M-P., Dukas, R., McClelland, G. B. (2010). Lifetime- and
caste-specific changes in flight metabolic rate and muiscle biochemistry
of honeybees, \emph{Apis mellifera} . \emph{J. Comp. Physiol. B.}
180(1): 45-55

Sim, C. \& Denlinger, D. L. (2013). Insulin signaling and the regulation
of insect diapause. \emph{Front. Physiol.} 4, 1--10

Soneson, C., Love, M. I., Robinson, M. D. (2015). Differential analysis
for RNA-seq: transcript-level estimates improve gene-level
inferences.\emph{F1000Research} , 4

Toth, A. L. \& Robinson, G. E. (2005). Worker nutrition and division of
labor in honey bees. \emph{Animal Behavior.} 69(2): 427-435

Toth, A. L., Kantarovich, S., Meisel, A. F., Robinson, G. E. (2005).
Nutritional status influences socially regulated foraging ontogeny in
honey bees. \emph{J. Exp. Bio.} 208: 4641-4649

Treanore, E. D., Kiner, J. M., Kerner, M. E., Amsalem, E. (2020). Shift
in worker physiology and gene expression pattern from reproductive to
diapause-like with colony age in the bumble bee \emph{Bombus impatiens}
.\emph{J. Exp. Biol.} 223, jeb218768

Vance, J. T., Williams, J. B., Elekonich, M. M., Roberts, S. P. (2009).
The effects of age and behavioral development on honey bee (\emph{Apis
mellifera} ) flight performance. \emph{J. Exp. Bio} . 212: 2604-2611

Wallberg, A. et al\emph{.} (2019). A hybrid de novo genome assembly of
the honeybee, Apis mellifera, with chromosome-length scaffolds.\emph{BMC
Genomics.} 20, 275

Whitfield, C. W., Cziko, A-M., Robinson, G. E. (2003). Gene expression
profiles in the brain predict behavior in individual honey
bees.\emph{Science} . 302(5643): 296-299

Winston, M. L. (1987). \emph{The Biology of the Honey Bee} . (Harvard
University Press)

Yamada, H. \& Yamamoto, M. T. (2011). Association between circadian
clock genes and diapause incidence in Drosophila triauraria. \emph{PLoS
ONE.} 6

Yocum, G. \& Rinehart, J. (2015). Key molecular processes of the
diapause to post-diapause quiescence transition in the alfalfa
leafcutting bee \emph{Megachile rotundata} identified by comparative
transcriptomic analysis. \emph{Physiol. Entomol.} 40, 103--112

Yocum, G. D. et al\emph{.} (2018). Environmental history impacts gene
expression during diapause development in the alfalfa leafcutting
bee,\emph{Megachile rotundata} . \emph{J. Exp. Biol.} 221, jeb173443

Yu-Xuan Lu et al. (2013). Polycomb repressive complex 2 (PRC2) protein
ESC regulates insect developmental timing by mediating H3K27me3 and
activating prothoracictropic hormone gene expression. \emph{J Biol
Chem.} 288(32): 23554-23564

Zhao, W. et al\emph{.} (2014). Comparison of RNA-Seq by Poly (A)
Capture, Ribosomal DNA

Depletion, and DNA Microarray for Expression Profiling. \emph{BMC
Genomics} . 15, 419

\textbf{Data Accessibility}

Transcriptomic data generated from this study have been deposited in the
NCBI Short Read Archive (SRA) database (BioProject accession
PRJNA601517). Reproducible source code is available at
https://sbresnahan.github.io/Winter-Honeybee-Transcriptome

\textbf{Author Contributions}

M.A.D., T.G. and C.M.G. designed the study. M.A.D. collected samples and
performed RNA extractions. S.T.B. performed transcriptomic and
statistical analyses. S.T.B. and C.M.G. prepared the initial drafts of
the manuscript and all authors contributed to the writing.

\textbf{Tables and Figures}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/Screen-Shot-2021-07-06-at-9.20.30-AM}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/Screen-Shot-2021-07-06-at-9.20.50-AM}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

(a) fat body sample correlations

\par\null\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig1a/fig1a}

\end{center}
\end{figure}

(b) flight muscle sample correlations

\par\null\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig1b/fig1b}

\end{center}
\end{figure}

\textbf{FIGURE 1} \textbf{(a)} fat body and \textbf{(b)} flight muscle
tissue sample relationships as determined by pairwise correlations and
hierarchical clustering. Sample relationships determined using all
tested genes (\emph{n =} 9,178).\textbf{Sample ID key} : Sample Group --
Tissue -- Replicate.

\par\null\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{FIGURE 2} \emph{k} -means clustering of \textbf{(a)} fat body
and\textbf{(b)} flight muscle tissue samples. Dotted lines indicate 95\%
confidence intervals.

\selectlanguage{english}
\FloatBarrier
\end{document}

