Differential gene expression underlying ovarian phenotype determination in honey bee, Apis mellifera L., caste development

Accepted Manuscript Differential gene expression underlying ovarian phenotype determination in honey bee, Apis mellifera L, caste development Denyse Cavalcante Lago, Fernanda Carvalho Humann, Angel Roberto Barchuk, Kuruvilla Joseph Abraham, Klaus Hartfelder PII:

S0965-1748(16)30137-0

DOI:

10.1016/j.ibmb.2016.10.001

Reference:

IB 2881

To appear in:

Insect Biochemistry and Molecular Biology

Received Date: 5 August 2016 Revised Date:

26 September 2016

Accepted Date: 3 October 2016

Please cite this article as: Lago, D.C., Humann, F.C., Barchuk, A.R., Abraham, K.J., Hartfelder, K., Differential gene expression underlying ovarian phenotype determination in honey bee, Apis mellifera L, caste development, Insect Biochemistry and Molecular Biology (2016), doi: 10.1016/j.ibmb.2016.10.001. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT 1

Differential gene expression underlying ovarian phenotype determination in

2

honey bee, Apis mellifera L, caste development

3

Denyse Cavalcante Lago*1, Fernanda Carvalho Humann*2,3, Angel Roberto Barchuk4,

4

Kuruvilla Joseph Abraham 5,6, Klaus Hartfelder2

5

1

6

de São Paulo, Avenida Bandeirantes 3900, 14049-900 Ribeirão Preto, SP, Brazil

7

2

8

de Medicina de Ribeirão Preto, Universidade de São Paulo, Avenida Bandeirantes

9

3900, 14049-900 Ribeirão Preto, SP, Brazil

RI PT

Departamento de Genética, Faculdade de Medicina de Ribeirão Preto, Universidade

SC

Departamento de Biologia Celular e Molecular e Bioagentes Patogênicos, Faculdade

10

3

11

Rua Estéfano D'avassi, 625, 15991-502 Matão, SP, Brazil

12

4

13

Biomédicas, Universidade Federal de Alfenas, Rua Gabriel Monteiro da Silva 700,

14

37130-000 Alfenas, MG, Brazil

15

5

16

Universidade de São Paulo, Avenida Bandeirantes 3900, 14049-900 Ribeirão Preto,

17

SP, Brazil

18

6

19

Preto, SP, Brazil

22

M AN U

TE D

Departamento de Puericultura e Pediatria Faculdade de Medicina de Ribeirão Preto,

Universidade Estácio-Uniseb, Rua Abrahão Issa Halach 980, 14096-160 Ribeirão

EP

21

Departamento de Biologia Celular e do Desenvolvimento, Instituto de Ciências

Short title: Gene expression in larval honey bee ovaries

AC C

20

Instituto Federal de Educação, Ciência e Tecnologia de São Paulo, Campus Matão,

23

* these authors contributed equally to the present study

24

Corresponding Author: Klaus Hartfelder, Departamento de Biologia Celular e

25

Molecular e Bioagentes Patogênicos, Faculdade de Medicina de Ribeirão Preto,

26

Universidade de São Paulo, Avenida Bandeirantes 3900, 14049-900 Ribeirão Preto,

27

SP, Brazil. E-mail: [email protected]

28 29

Co-author e-mail addresses: [email protected] (DCL); fchumann@gmail,com

30

(FCH); [email protected] (ARB); [email protected] (KJA)

31

1

ACCEPTED MANUSCRIPT ABSTRACT

33

Adult honey bee queens and workers drastically differ in ovary size. This adult ovary

34

phenotype difference becomes established during the final larval instars, when

35

massive programmed cell death leads to the degeneration of 95-99% of the ovariole

36

anlagen in workers. The higher juvenile hormone (JH) levels in queen larvae protect

37

the ovaries against such degeneration. To gain insights into the molecular architecture

38

underlying this divergence critical for adult caste fate and worker sterility, we

39

performed a microarray analysis on fourth and early fifth instar queen and worker

40

ovaries. For the fourth instar we found nine differentially expressed genes (DEGs)

41

with log2FC > 1.0, but this number increased to 56 in early fifth-instar ovaries. We

42

selected 15 DEGs for quantitative PCR (RT-qPCR) analysis. Nine differed

43

significantly by the variables caste and/or development. Interestingly, genes with

44

enzyme functions were higher expressed in workers, while those related to

45

transcription and signaling had higher transcript levels in queens. For the RT-qPCR

46

confirmed genes we analyzed their response to JH. This revealed a significant up-

47

regulation for two genes, a short chain dehydrogenase reductase (sdr) and a heat

48

shock protein 90 (hsp90). Five other genes, including hsp60 and hexamerin 70b

49

(hex70b), were significantly down-regulated by JH. The sdr gene had previously

50

come up as differentially expressed in other transcriptome analyses on honey bee

51

larvae and heat shock proteins are frequently involved in insect hormone responses,

52

this making them interesting candidates for further functional assays.

AC C

EP

TE D

M AN U

SC

RI PT

32

53 54

Keywords: honeybee; programmed cell death; ovary transcriptome; juvenile

55

hormone; short chain dehydrogenase reductase; hexamerin

56

2

ACCEPTED MANUSCRIPT 57

1. Introduction Social insects, which can represent up to 30% of the animal biomass in

59

terrestrial ecosystems, are ecological success models (Fittkau and Kilinge, 1973;

60

Hölldobler and Wilson, 1990). This success story is based on reproductive division of

61

labor within a caste system, which, in the Hymenoptera, comprises a highly fertile

62

queen and a facultatively or permanently sterile worker caste. While these castes are

63

well-studied in terms of their ecology and behavior, understanding the molecular and

64

cellular processes underlying the development of the queen and worker phenotypes is

65

still a challenge.

M AN U

SC

RI PT

58

In the honey bee, Apis mellifera, which is a model system for social

67

organization even beyond the insects, a queen typically has 150-200 ovarioles in each

68

of her ovaries and is capable of laying over 1,000 eggs per day, whereas the worker

69

ovary typically consists of only 2-12 of these serial reproductive system units

70

(Martins and Serrão, 2004; Rozen, 1986). Furthermore, honey bee workers typically

71

do not lay eggs in the presence of a queen. Despite the wide gap between queen and

72

worker ovariole numbers, the actual number of ovarioles is, however, quite variable

73

within each of the two castes. In workers, such variance was shown to be associated

74

with the paternal genotype (Makert et al., 2006), and several quantitative trait loci

75

affecting worker ovariole number have been revealed through backcrosses of African

76

and European honey bee stocks (Linksvayer et al., 2009; Linksvayer et al., 2011). In

77

contrast, the gap in ovary size (ovariole number) between queens and workers is

78

primarily linked to the nutritional conditions they experienced during larval

79

development (Dedej et al., 1998). Nonetheless, the binomial ovary size distribution

80

that separates naturally reared queens from natural workers essentially transforms into

81

a continuum when larvae are reared in vitro in the laboratory (Leimar et al., 2012).

AC C

EP

TE D

66

3

ACCEPTED MANUSCRIPT The question then is, what are the factors that generate this binomial gap seen under

83

the “natural experiment” conditions. Initially triggered by the diets that the larvae

84

receive from nurse bees, morphogenetic hormones and nutrient sensing signaling

85

pathways become then critical factors for caste phenotype decisions in honey bees

86

(for a recent review see Hartfelder et al., 2015).

RI PT

82

While this clearly is a question about proximate ontogenetic mechanisms, the

88

exaggerated ovary size seen in queens of the genus Apis is also of interest in terms of

89

evolutionary developmental biology. Even though ovariole numbers are highly

90

variable across and within insect orders (see Büning, 1994 for a comprehensive

91

compilation of insect ovary character states), three ovarioles per ovary appear to be

92

the basal state for bees (Iwata, 1955; Martins and Serrão, 2004; Rozen, 1986). In the

93

four tribes that form a monophyletic clade within the subfamily Apinae, and which

94

include the Apini (Cardinal et al., 2010; Michener, 2007), four ovarioles per ovary are

95

typical for orchid bees (Euglossini) and bumble bees (Bombini). Ovariole number

96

also appears to be fixed to four in the highly eusocial stingless bees (Meliponini),

97

although some variation in ovariole number has been observed in queens (Cruz-

98

Landim et al., 1998). Thus, the genus Apis, with its wide gap in ovariole numbers

99

between the two castes, is an outstanding case, not only within this clade but even

100

across all hymenopteran orders. It is only championed by queens of army ants, where

101

around 250 ovarioles per ovary are reported for Eciton schmitti queens (Wheeler,

102

1910), and this is even further surpassed by the extreme case of queens of the African

103

driver ant genus Dorylus, which can have up to 15,000 ovarioles, making them

104

capable of laying 1 to 2 million eggs per month (Hölldobler and Wilson, 1990).

AC C

EP

TE D

M AN U

SC

87

105

Exaggerated morphologies are by no means limited to the caste phenotypes in

106

social insects, but have evolved in many insect orders that exhibit marked genetic or

4

ACCEPTED MANUSCRIPT facultative polyphenisms within or among sexes, and there appear to be common

108

themes in terms of developmental regulation via endocrine system functions

109

associated with metamorphosis (Hartfelder and Emlen, 2012). In the honey bee, one

110

of the major effectors in caste fate decision is the larval juvenile hormone (JH) titer,

111

which is low in worker and high in queen larvae up to the beginning of the fifth instar

112

(Rachinsky et al., 1990; Rembold, 1987). The elevated JH level in the hemolymph of

113

queen larvae is crucial for preventing the induction of massive programmed cell death

114

in their ovaries, while 95-99% of the ovariole primordia are degraded in fifth-instar

115

worker larvae (Schmidt Capella and Hartfelder, 1998).

SC

RI PT

107

With the aim of revealing molecular mechanisms and pathways underlying the

117

caste-specific developmental trajectories of the honey bee ovary we performed a

118

microarray analysis on ovaries of fourth and early fifth-instar queen and worker

119

larvae, i.e. the period during which we expected transcriptional decisions critical for

120

ovarian phenotype determination. We found a clear increase in the number of

121

differentially expressed genes (DEGs) soon after the larvae had entered the final fifth

122

instar. A set of 15 genes indicated as differentially expressed in the microarray

123

analyses was further analyzed by real-time PCR and assayed with respect to their

124

response to juvenile hormone.

TE D

EP

AC C

125

M AN U

116

126

2. Material and methods

127

2.1. Honey bee larval ovary sampling

128

Fourth (L4) and feeding phase fifth (L5F1-L5F3) instar queen and worker

129

larvae were obtained from hives of Africanized hybrid Apis mellifera stocks kept in

130

the Experimental Apiary of the University of São Paulo in Ribeirão Preto. The larvae

5

ACCEPTED MANUSCRIPT were staged according to established criteria (Michelette and Soares, 1993; Rachinsky

132

et al., 1990). Their ovaries were dissected in sterile honey bee saline (Rachinsky et al.,

133

1990) and immediately transferred to TRIzol reagent (Life Technologies, Carlsbad,

134

CA, USA) for storage at -80 °C. RNA extraction then followed the manufacturer’s

135

protocol (Life Technologies), with glycogen as an additive to improve pellet

136

formation.

RI PT

131

For the microarray samples, 40 ovary pairs each were dissected from fourth

138

(L4) and early fifth (L5F1/L5F2) queen and worker instar larvae (total number of

139

samples, N = 4). For the RT-qPCR analyses of DEGs and the respective JH

140

responses, 10 ovary pairs were dissected for each biological sample (total number of

141

samples, N = 33).

142

M AN U

SC

137

2.2. Microarray analysis

144

2.2.1. Sample preparation and data acquisition

TE D

143

RNA extracted from the larval ovaries was further purified using the RNA

146

Cleanup kit (RNeasy Mini Kit, Qiagen, Hilden, Germany), and 1 µg of each total

147

RNA was then used for amplification with the Amino Allyl MessageAmp™ II aRNA

148

Amplification Kit (Ambion, Foster City, CA, USA), following the manufacturer’s

149

instructions. For probe preparation in a dye-swap hybridization design, 20 µg of each

150

amplified amino allyl RNA was subsequently labeled with Cy3 or Cy5 dyes (CY Dye

151

Post-Labeling Reactive Dye Pack, GE Healthcare Life Sciences, Pittsburgh, PA,

152

USA).

AC C

EP

145

153

Mutual sets of Cy3/Cy5 labeled probes were hybridized to honey bee whole-

154

genome oligonucleotide microarrays acquired from the Functional Genomics Unit of 6

ACCEPTED MANUSCRIPT 155

the W.M. Keck Center at the University of Illinois at Urbana-Champaign. The slide

156

design

157

microarray hybridization process followed the protocol described by Moda et al.

158

(2013) established for honey bees.

159

2.2.2. Microarray – data analysis

available

at:

http://www.biotech.uiuc.edu/functionalgenomics.

The

RI PT

is

The scan results of the slides read in an Axon Genepix 4000B scanner

161

(Molecular Devices, Sunnyvale, CA, USA) with GenePix® Pro 6.0 software (Agilent

162

Technologies, Santa Clara, CA) were analyzed using the Limma package of

163

Bioconductor software (Gentleman et al., 2004) implemented in R (R Development

164

Core Team, 2011). The Limma package was used for normalization between and

165

within arrays for each comparison. DEGs were then identified by fitting a Linear

166

Model, followed by empirical Bayes methods for determining variances.

M AN U

SC

160

The data analysis strategy followed the protocol described by Moda et al.

168

(2013). Briefly, background correction was done by adding a positive constant (offset

169

correction = 50) to the background intensities. This restrains spurious variation in log

170

ratios. For correcting within-assay dye intensity variation and spatial effects, a “print-

171

tip loess” normalization was used, and single channel normalization was done to

172

facilitate comparisons across arrays. Subsequently, the log2 ratio for queen versus

173

worker sample intensity was determined for each probe in each array. After applying

174

Bayesian statistics to shrink estimated standard errors, a moderated t-statistic test was

175

run for each gene. This analysis permitted to score DEGs of queen versus worker

176

samples by their log2 fold change (log2 FC) values and B probability. The data were

177

deposited in the Gene Expression Omnibus database (GEO, at the NCBI

178

database) under accession number GSE81118.

AC C

EP

TE D

167

7

ACCEPTED MANUSCRIPT Next, to check for eventual false positives that might arise from multiple

180

testing, the Benjamini-Hochberg method was employed, with an adjusted P-value

181

threshold of 0.05, but no threshold was imposed on Fold Change. With this criterion

182

for identifying DEGs, 30 top ranking genes from the L5 ovary analysis were selected

183

for a heat map and cluster analysis using the gplots2 package in R (R Development

184

Core Team, 2011).

RI PT

179

185

187

SC

186

2.3. Real-time quantitative PCR analysis of differentially expressed genes (DEGs) For an in-depth analysis on differential gene expression in the ovaries of queen

189

and worker larvae we selected 15 genes from the list indicated by the microarray

190

analysis. The primary criterion for inclusion was a log2FC value ≥1.5. In the case of

191

genes annotated with potential functions in cell signaling pathways, a log2FC value

192

>1.0<1.5 was admitted. These 15 genes and their respective primers are listed in

193

Table 1. Primers

were

TE D

194

M AN U

188

designed

using

the

Primer-BLAST

tool

at

NCBI

(http://www.ncbi.nlm.nih.gov/tools/primer-blast/). Primer conditions and amplicon

196

specificity were first checked by PCR in a gradient thermocycler (PTC-200, MJ

197

Research, Waltham, MA, USA). The amplification products were electrophoretically

198

separated in agarose gels for product size confirmation. For confirmation of PCR

199

product identity these were purified using the Wizard SV Gel and PCR Clean-Up

200

System (Promega, Madison Wi, USA) and either directly sequenced, or cloned by

201

insertion into a plasmidial vector (pGEM-T Easy Vector Systems, Promega) and

202

transformation of E. coli DH5α bacteria. Sequencing was done using the Big Dye

203

Terminator Cycle Sequencing Ready Reaction (Applied Biosystems, Foster City, CA,

AC C

EP

195

8

ACCEPTED MANUSCRIPT 204

USA) and M13 forward primer on an ABI-PRISM 3100 (Applied Biosystems)

205

automated gene analyzer. The obtained sequences were compared with those of the

206

annotated honey bee genes of interest using the BLASTN algorithm. TRIzol-extracted RNA from dissected ovaries was treated with 0.1 U DNaseI

208

(Invitrogen) and RNA quality and quantity were assessed spectrophotometrically.

209

First strand cDNA was produced using a SuperScript II (Invitrogen) protocol at 42°C

210

for 50 min and 70°C for 15 min. Quantitative RT-PCR (RT-qPCR) analyses were set

211

up using 1.5 µL cDNA (diluted 1:10), 7.5 µL of Power SYBR® PCR Green Master

212

Mix (Life Technologies), 0.5 µL of each forward and reverse primer (10 pmol/µL)

213

and 5 µL of deiononized water (MilliQ, Millipore, Billerica, MA, USA), totalizing a

214

final volume of 15 µL. Reactions were run in a Real-Time PCR StepOne Plus system

215

(Applied Biosystems) under the following conditions: 50 °C for 2 min, 95 °C for 10

216

min, 40 cycles of 95 °C for 15 s and 60 °C for 1 min, followed by dissociation curve

217

analysis (95 °C for 15 s, 60 °C for 1 min and 95 °C for 15 s). For normalization of

218

the biological samples we performed assays using the honey bee rp49 gene (BeeBase

219

code in v3.2, GB47227), which encodes a ribosomal protein (RP49, also named

220

RPL32), and a gene encoding a cytoplasmic actin (BeeBase code in v3.2, GB44311).

221

Their suitability as endogenous control genes in A. mellifera RT-qPCR assays has

222

previously been established (Lourenço et al., 2008). Relative expression was

223

calculated following the method of Livak and Schmittgen (2001).

SC

M AN U

TE D

EP

AC C

224

RI PT

207

To assess primer efficiencies for RT-qPCR, cDNA from six samples of queen

225

and worker larval ovaries were pooled and used to prepare a 10x-step dilution series.

226

The Ct results for these series were used to calculate slope, r2, dissociation

227

temperature and efficiency (%) for each primer set. Only primer sets with efficiencies

228

between 80% and 105% were used in the subsequent RT-qPCR assays, which were

9

ACCEPTED MANUSCRIPT 229

all done as triplicate runs for each biological sample, and three biological samples

230

were analyzed for each developmental stage and caste.

231

2.4. Juvenile hormone treatment effect on gene expression in larval worker ovaries

RI PT

232

This experiment was done to assess whether JH may directly affect the

234

expression of certain genes for which caste differences in transcript levels had been

235

confirmed by the previous RT-qPCR assays. For the treatment we retrieved a brood

236

frame containing fourth-instar larvae and applied a 10 µg dose of juvenile hormone

237

III (Fluka, Buchs, Switzerland) dissolved in acetone (10 µg/µL) directly on their

238

cuticle. Control larvae were either left untreated or treated with 1 µL of acetone

239

(solvent control). The frame was returned to the colony of origin, and 6 h after the

240

treatment the larvae were removed and had their ovaries dissected for cDNA

241

preparation.

TE D

M AN U

SC

233

In the RT-qPCR analyses we assayed the transcript levels of the following

243

genes (BeeBase codes are those for version Official Gene set version 3.2): apoLp-III

244

(GB55452), gpdh (GB50902), hex70b (GB51697), hsp60 (GB54372) and hsp90

245

(GB40976), mapk-3 (GB41845), oclp-1 (GB55030), 15-pgdh (GB46366) and sdr

246

(GB54419), all of which had been confirmed as differentially expressed in the prior

247

time course analysis. So as to confirm the efficiency of the JH treatment, we also

248

analyzed the transcript levels of the honey bee krüppel homolog-1 encoding gene (kr-

249

h1; GB45427; XM_001123084), which is an immediate response gene in the JH

250

response cascade (Belles and Santos, 2014; Jindra et al., 2013).

AC C

EP

242

251

10

ACCEPTED MANUSCRIPT 252

2.5. Statistical analysis Transcript levels for the 15 genes selected from the microarray results were

254

compared with respect to caste and developmental stage by means of two-way

255

ANOVA, followed by a Sidak multiple comparisons post hoc test. The analysis of JH

256

treatment effects was done using a Brown-Forsythe test for assessing equality of

257

variance and one-way ANOVA followed by a Holm-Sidak multiple comparisons test.

258

The statistical analyses were done using GraphPad Prism v. 5.0 software.

SC

RI PT

253

M AN U

259

260

3. Results

261

3.1. Microarray analysis of differential gene expression in the honey bee larval

262

ovary

For this analysis we focused on the fourth and early fifth larval instar, as it is

264

during this period that the ovary phenotypes are strongly affected by the larval rearing

265

environment (Dedej et al., 1998). These two stages also comprise the period just prior

266

to or coinciding with the beginning of massive programmed cell death in the larval

267

worker ovary (Schmidt Capella and Hartfelder, 1998; Hartfelder and Steinbrück,

268

1997; Reginato and Cruz-Landim, 2001). In the dye-swap design, the RNA extracted

269

from the ovaries generated two microarray sets each, one contrasting fourth and one

270

early fifth-instar ovaries for each caste. Genes were considered as differentially

271

expressed when a log2-fold change (log2 FC) >1.0 was attained. Table 2 lists the genes

272

that were more highly expressed in larval queen ovaries, and Table 3 shows those that

273

were higher expressed in worker ovaries

AC C

EP

TE D

263

11

ACCEPTED MANUSCRIPT In general terms it was possible to see that only very few genes were classified

275

as differentially expressed (log2FC ≥ 1) in the L4 ovary microarrays, two as higher

276

expressed in queen (Table 2) and seven in worker ovaries (Table 3). This changed

277

drastically once the larvae had entered the fifth instar, where 32 microarray spots gave

278

higher expression signals for queens (Table 2) and 24 for workers (Table 3). These

279

molecular signatures can be interpreted as an incremental increase in gene expression

280

differences from a still bipotential ovary phenotype in fourth-instar larvae to a gradual

281

fixation of ovarian caste fate in the early fifth instar.

SC

RI PT

274

To further investigate this transition in gene expression patterns from the

283

fourth to the fifth larval instar and to avoid eventual pitfalls due to false positives we

284

ran a second statistical analysis protocol for the microarray data using a Benjamini-

285

Hochberg correction for multiple testing with a threshold of p<0.05, but without

286

setting a threshold for fold change. This generated a second gene list of 30 genes for

287

the L5 ovary microarrays that included 24 of the ones listed in Tables 2 and 3. The

288

differential expression results for these top 30 genes are shown in a heat map

289

representation (Figure 1A), where the fourth (L4) and fifth instar (L5) ovaries were

290

represented as duplicated by their respective dye-swap reads. This illustrates that the

291

fifth-instar worker ovaries are divergent in their gene expression profile not only from

292

fifth-instar queen but also from fourth-instar ovaries of both castes, thus underlining

293

the importance of the entry into the last fifth instar as a critical time point for caste

294

fate determination of the worker ovary phenotype.

AC C

EP

TE D

M AN U

282

295

When comparing the gene lists in Tables 2 and 3 in terms of general

296

characteristics, the two lists turned out to be rather similar. Among the DEGs that are

297

more highly expressed in queen ovaries (Table 2), 19 are protein-coding genes of

298

known or structure-predicted function (Apidaecin 14 and Apidaecin 22 were counted

12

ACCEPTED MANUSCRIPT only once each), and seven ESTs represent sequences of unknown function. Among

300

the DEGs that are more highly expressed in worker ovaries (Table 3), 16 are protein-

301

coding genes of known or structure-predicted function (hexamerin 70c and serine

302

protease 17 were only counted once each), four genes represent hypothetical proteins,

303

and six are ESTs. Yet, when comparing the DEGs in terms of functional categories

304

(Figure 1B), those listed as higher expressed in queen ovaries (Table 2) comprised

305

more genes related to expression (transcription, translation) and regulation, as well as

306

to defense and immune response. In contrast, among the genes indicated as higher

307

expressed in worker ovaries (Table 3) there were more annotated as related to

308

transport functions and enzyme inhibition (but note, no statistical test on GO

309

categories was done due to the low number of DEGs).

M AN U

SC

RI PT

299

310

3.2. Transcript levels of DEGs during larval ovary development

TE D

311

So as to strengthen the conclusions on differential gene expression patterns

313

revealed through the microarray results we selected a set of 18 genes from those listed

314

in Tables 2 and 3 for an in-depth time course study of their expression. The temporal

315

window covered the fourth larval (L4) instar and the three feeding substages of the

316

fifth larval instar (L5F1-L5F3), these representing the stages during which the larvae

317

are still being fed by nurse bees and grow in size. Furthermore, it was in the last of

318

these stages (L5F3) that the onset of massive programmed cell death had been

319

revealed by ultrastructure analysis (Hartfelder and Steinbrück, 1997) and TUNEL

320

labeling (Schmidt Capella and Hartfelder, 1998). For 15 of these 18 genes we could

321

design primers of adequate efficiency for relative quantification of expression by real-

322

time PCR (Table 1), and, by two-way ANOVA analysis, six of these were found to

AC C

EP

312

13

ACCEPTED MANUSCRIPT 323

differ significantly in their transcript levels with respect to caste and/or developmental

324

stage (Figure 2). Functionally, these genes could be grouped into the following functional

326

categories: enzymes, represented by four genes, factors involved in transcription and

327

cellular signaling (six genes), storage functions (two genes), defense (two genes), and

328

a gene with an associated hypothetical function. Among the enzyme encoding genes

329

(Figure 2A), three had consistently elevated mean transcript levels in worker ovaries,

330

these being a short-chain reductase hydrogenase (sdr), a 15-hydroxy prostaglandin

331

dehydrogenase (15-pgdh) and a retinoid-inducible carboxypeptidase-like enzyme

332

(acpep1). The fourth gene in the enzyme category, a glycerol-3-phosphate

333

dehydrogenase (gpdh), was also more highly expressed in the L5F1 stage, but then the

334

pattern became inverted, and it was significantly overexpressed about 24 h later in

335

ovaries of the L5F3 stage. Nonetheless, statistical significance for differential

336

expression could only be confirmed for the sdr gene in the L5F1 stage, and the gpdh

337

gene in the L5F3 stage (for details on statistical analysis results see Supplementary

338

material Table S1).

TE D

M AN U

SC

RI PT

325

Among the six genes related to transcription and cellular signaling functions

340

(Figure 2B), two, the heat shock protein 60 (hsp60), heat shock protein 90 (hsp90)

341

genes, were significantly higher expressed in queen larval ovaries. No significant

342

differences were found for the mitogen-activated protein kinase 3 (mapk-3), helicase

343

25E (hel25e), odorant binding protein 14 (obp14) and the potential ecdysone-induced

344

protein 93 (e93) gene, which is also known as mblk-1 in the honey bee.

AC C

EP

339

345

With respect to storage function-related proteins (Figure 2C), the

346

apolipoprotein III (apoLp-III) gene showed increasingly higher expression levels in

347

worker ovaries, whereas the hexamerin 70b (hex70b) had higher mean transcript

14

ACCEPTED MANUSCRIPT levels in L5F1 queen ovaries. Among the two defense functions-related proteins, the

349

hypothetical inhibitor function cysteine knot peptide gene (oclp-1) had significantly

350

higher expression levels in L5F1- and L5F3-stage worker ovaries (Figure 2D). The

351

differential expression for apidaecin 14 (apid14) was not confirmed by the RT-qPCR

352

assays. For the gene encoding a hypothetical protein (GB43690) no significant

353

expression difference was found with respect to caste and larval stage.

RI PT

348

When looking from a general perspective we noted that the directionality of

355

differential expression between the microarray and the RT-qPCR data appeared

356

consistent for several but not all of the genes. This may, in part, be explained by the

357

facts that: 1) the microarrays were done with only one biological sample each for the

358

fourth and early fifth instar of each caste, even though these were composed of 40

359

ovaries each; 2) different biological samples were used in the two assay types; and 3)

360

ovaries from the three substages of the fifth-instar feeding phase had been pooled in

361

the microarrays. Hence, for an overall comparison of the two datasets we first

362

calculated mean log2FC values for the three fifth-instar substages for the genes

363

analyzed by RT-qPCR and then plotted these against the respective microarray

364

log2FC values (Figure 1S, Supplementary material). The linear regression with an

365

r2=0.3492 was significant (F1,12=6.483, p=0.261) and so was the subsequent Spearman

366

rank correlation analysis (p=0.0209). This confirms that there is good overall

367

comparability across the two datasets.

369

M AN U

TE D

EP

AC C

368

SC

354

3.3. Juvenile hormone effect on transcript levels

370

Since juvenile hormone plays an important role in defining the adult ovary

371

phenotype by inhibiting programmed cell death during larval development (Schmidt

15

ACCEPTED MANUSCRIPT Capella and Hartfelder, 1998) we next checked whether this hormone affects the

373

expression of the genes confirmed as differentially expressed by real-time PCR. A JH

374

dose as low as 1 µg is known to be effective in inducing queen-like morphological

375

characters when topically applied to honey bee worker larvae (Dietz et al., 1979;

376

Rembold et al., 1974), and treatment of larvae with a 10 µg dose of JH III has been

377

shown to inhibit programmed cell death in worker ovaries (Schmidt Capella and

378

Hartfelder, 1998). Since we chose a very narrow time window of only six hours to

379

reveal direct gene regulatory effects of this hormone in larval ovaries, we decided to

380

use this higher dose. As a positive control for the JH response, we assayed the

381

expression of the immediate response gene kr-h1 (Belles and Santos, 2014; Jindra et

382

al., 2013).

M AN U

SC

RI PT

372

As expected, kr-h1 expression in ovaries of early fifth-instar worker larvae

384

was significantly up-regulated by the topical JH application (Figure 3A; F2,8=46.09;

385

p<0.0001) in comparison to untreated control larvae. There was also a clear acetone

386

solvent effect, but this went in the opposite direction, thus further putting in evidence

387

the strong effect of JH on the induction of this immediate response gene. Two genes,

388

sdr and hsp90 were significantly up-regulated by the JH application, but interestingly,

389

their expression was not affected by the solvent (sdr F2,6=2.515, p=0.1610; hsp90

390

F2,8=5.898, p=0.0267) .

EP

AC C

391

TE D

383

An opposite effect of JH was seen for 15-pgdh, oclp-1, apoLp-III, hex70b and

392

hsp60. All of these genes were significantly down-regulated by JH in comparison to

393

the untreated and solvent-treated controls (Figure 3B; 15-pgdh F2,8=7.021, p=0.0174;

394

oclp-1 F2,8=6.938 p=0.0179; apoLp-III F2,8=6.661 p=0.0198; hex70b F2,8=8.123

395

p=0.0119; hsp60 F2,8=6.661 p=0.0198). Furthermore, the expression of all these five

396

genes appeared to be affected by the solvent treatment, but not strongly enough to be

16

ACCEPTED MANUSCRIPT 397

significant. The other two genes included in this analysis, mapk-3 and gpdh, were not

398

affected, neither by JH nor the solvent (data not shown).

399

4. Discussion

RI PT

400

Being a model organism for genomic studies on insect sociality, gene

402

expression analyses in the honey bee have traditionally set their focus on three

403

aspects, division of labor among adult workers, reproduction and reproductive

404

conflict among adult bees, and caste development in the larval stages. This is

405

especially the case for large-scale transcriptome studies. The first comprehensive one,

406

based on EST libraries and done even before the whole genome had been sequenced

407

and annotated, was a microarray comparison of the brain transcriptomes of young

408

(nurse) and old (forager) worker bees (Whitfield et al., 2003). Head transcriptomes of

409

nurses and foragers were subsequently analyzed by RNA-seq (Liu et al., 2013). The

410

question of why adult workers normally do not reproduce in the presence of a queen,

411

i.e. facultative worker sterility, was assessed by microarray comparisons of abdominal

412

and brain RNA of workers exhibiting different degrees of ovary activation

413

(Thompson et al., 2006), and subsequently through comparative ovary and abdomen

414

transcriptome analysis of anarchist versus wild-type worker bees (Kucharski et al.,

415

2008). More recently, an RNA-seq analysis on queen and worker ovaries of different

416

activation states revealed high numbers of differentially expressed genes, especially

417

for inactive versus active ovaries of workers (Niu et al., 2014).

AC C

EP

TE D

M AN U

SC

401

418

For caste development, the first microarray comparison was also based on an

419

EST set, showing that genes related to developmental regulation were higher

420

expressed in worker larvae, while those related to incremental growth were more

17

ACCEPTED MANUSCRIPT expressed in queens (Barchuk et al., 2007). Subsequently, the development of a

422

genome-based oligonucleotide-spotted microarray platform for A. mellifera (W.M.

423

Keck Center, University of Illinois) enlarged and standardized the microarray

424

analyses for this species. This microarray platform was directly compared to Illumina

425

RNA-seq data in a study on early divergence of gene expression for honey bee queen

426

and worker larvae (Cameron et al., 2013). The study revealed distinct patterns of

427

DEGs for the early, mid and late larval stages and identified major classes associated

428

with caste, such as cytochrome P450 and spliceosome factors encoding genes. Yet,

429

despite being highly informative and potential material for future gene interaction and

430

regulatory network analysis, the analysis was done on whole body RNA extracts, and

431

knowingly, in the larval stages, different tissues have different developmental

432

dynamics. For instance, the major body tissue, the larval fat body, is the center of the

433

intermediary and energy metabolism and will be completely destroyed and

434

reorganized during metamorphosis, while others, such as the leg and wing imaginal

435

discs are composed of mitotic cells that have low metabolic rates and undergo

436

differentiation only during metamorphosis.

TE D

M AN U

SC

RI PT

421

Like the imaginal discs that form thoracic appendages, the ovary is also an

438

organ that will attain its functional capacity only in the adult stage. Furthermore, as it

439

is only a small compartment in the larval body, its contribution to DEGs in whole

440

body analyses of larvae is likely minimal, compared to its fundamental importance for

441

the adult caste phenotype. In a previous study, where we used a Representational

442

Difference Analysis to enrich for DEGs of queen and worker ovaries from final

443

feeding stage larvae (L5F3), we sequenced and identified over 70 ESTs (Humann and

444

Hartfelder, 2011). Among these, two annotated as potential long noncoding RNAs

445

(Humann et al., 2013) were further analyzed by RT-qPCR and shown to map to a

AC C

EP

437

18

ACCEPTED MANUSCRIPT QTL for variation in worker ovariole number (Linksvayer et al., 2011). With the

447

current microarray-based transcriptome analysis we expected to gain a more ample

448

picture on DEGs for queen and worker ovaries. Furthermore, we extended the

449

analysis to developmental stages that precede the onset of massive programmed cell

450

death, which starts to become manifest in L5F3 worker ovaries (Schmidt Capella and

451

Hartfelder, 1998).

452 4.1. DEGs in queen versus worker honey bee larval ovaries

SC

453

RI PT

446

The microarray analysis put in evidence that the molt from the fourth to the

455

fifth instar is associated with a major transition in the gene expression profile of the

456

larval ovary. The number of DEGs (log2FC ≥ 1, in Tables 2 and 3) showed an

457

increase from nine to 47 (not counting duplicate representations of genes). This

458

transition is also well documented in the heatmap (Figure 1A) generated from the

459

microarray data. Interestingly, the heatmap representation revealed that the fifth-instar

460

worker ovaries are the most distinct ones in terms of gene expression, both from

461

queen fifth-instar ovaries, as well as from the fourth-instar ones, a fact that could be

462

related to the onset of the cell death program (Hartfelder and Steinbrück, 1997;

463

Reginato and Cruz-Landim, 2001).

TE D

EP

AC C

464

M AN U

454

For more detailed information on the expression of specific genes we next

465

performed RT-qPCR assays on 15 genes selected based on differential expression

466

levels or function. Differential expression and concordance in directionality (higher in

467

queen or worker) was confirmed for seven genes, hsp60 and hsp90, gpdh, acpep1,

468

mapk-3, oc1p-1 and apo-Lp-III. Differential expression, but with inverted

469

directionality was seen for hex70b. For three genes, differential expression indicated

470

from the microarrays was not confirmed (e93, hel25e and GB43690). For two genes,

19

ACCEPTED MANUSCRIPT sdr and obp14, which had appeared in both lists (Tables 2 and 3) as higher expressed,

472

but with different timing, no significant differential expression was seen for obp14 in

473

the RT-qPCR assays. But for sdr we saw a significantly higher expression in L5F1

474

worker ovaries, contrary to what the microarrays had indicated. Nonetheless, the

475

correlation analysis (Supplementary Material Figure 1S) done after correcting for

476

differences in sampling characteristics showed a significant coincidence for the

477

microarray with the RT-qPCR data.

RI PT

471

The fact there were no specific cell death-related genes in the DEG lists may

479

have two possible explanation. One is that the larval stages we looked at were those

480

prior to the most extensive phase of programmed cell death previously revealed by

481

ultrastructure and TUNEL analysis (Hartfelder and Steinbrück, 1997; Schmidt

482

Capella and Hartfelder, 1998). For instance, the expression of specific cell death

483

genes in honey bee ovaries has been investigated by Dallacqua and Bitondi (2014)

484

showing that an antiapoptotic gene (Ambuffy) and an apoptosis promoting gene

485

(Amark) become strongly expressed only after the larvae had stopped feeding and

486

entered metamorphosis. A second explanation is that the actual cell death executing

487

factors, the caspases, are all sequentially activated by proteolytic cleavage and, hence,

488

are not necessarily differentially expressed. Nonetheless, some of the DEGs related to

489

signaling pathways, such as mapk-3, hsp60 and hsp90 could represent early cell

490

physiology responses to nutrition and/or stress. Another DEG, which we did not study

491

further here, encodes a Tudor staphylococcal nuclease. This is an ancient family of

492

proteins involved in multiple cell functions, including stress responses, and in plants it

493

substitutes caspases (Gutierrez-Beltran et al., 2016). Nonetheless, all these genes play

494

multiple roles in cell cycle regulation and transcriptional control and are not

495

specifically associated with programmed cell death.

AC C

EP

TE D

M AN U

SC

478

20

ACCEPTED MANUSCRIPT 496 497

4.2. Juvenile hormone response of DEGs in worker larval ovaries Since JH is an important factor for caste-specific ovary development we asked

499

how the RT-qPCR confirmed DEGs might respond to this hormone. In the L4 and

500

L5F1 stages, the JH titer of queen larvae is approximately 10 times higher than in

501

worker larvae, but then the JH levels in queen larvae drop to levels similar to those in

502

workers as the larvae finish feeding and prepare for metamorphosis (Rachinsky et al.,

503

1990). With this in mind, we topically applied a known effective dose of JH III to

504

L5F1 worker larvae, and these were then sampled six hours later for RNA extraction

505

from the ovaries. We chose this short time window between hormone application and

506

sampling because we were primarily interested in relatively direct transcriptional

507

effects on the genes of interest. The suitability of this time window and also the

508

efficiency of the treatment were confirmed by quantification of the transcript levels

509

for the JH immediate response gene kr-h1. Of the 10 genes selected for this

510

experiment, eight showed a significant alteration in their transcript levels in response

511

to JH treatment. Three of these, hsp90, oclp1 and 15-pdgh, responded in the expected

512

directionality: hsp90, which is higher expressed in queens was up-regulated by JH in

513

worker ovaries; in contrast, oclp-1 and 15-pdgh, both being higher expressed in

514

worker ovaries, were down-regulated by the elevated JH levels. The other three genes,

515

sdr, hex70b and hsp60, showed a response opposite to expectation: sdr, which is

516

higher expressed in worker L5F1 ovaries was up-regulated by JH, and hex70b and

517

hsp60, which are higher expressed in L5F1 queen ovaries, were down-regulated by

518

JH. When taking together the results from the three analyses, microarray, RT-qPCR

519

and JH response, it is possible to say that hsp90 could be a molecular marker for

AC C

EP

TE D

M AN U

SC

RI PT

498

21

ACCEPTED MANUSCRIPT 520

queen-like ovary development, while 15-pgdh, oclp-1 and apoLp-III could be

521

candidate markers for worker-like ovary development.

522 523

4.3. Candidate marker genes for caste-specific ovary development Different from the candidate marker genes for worker ovary development, the

525

queen marker hsp90 consistently appeared as 4-5 times more expressed in queen

526

ovaries and was also seen to be strongly up-regulated by JH. The Hsp83 protein

527

homolog in the fruit fly has been shown to facilitate the nuclear import of the JH

528

receptor Methoprene-tolerant/Germ-cell-expressed, thus mediating the direct cellular

529

response to circulating JH levels (He et al., 2014). Furthermore, insect Hsp90 proteins

530

are important mediators of the crosstalk between the JH and ecdysteroid signaling

531

pathways (Cai et al., 2014; Liu et al., 2013). By playing such an important role in the

532

endocrine control mechanisms of insect growth and metamorphosis, this gene should

533

clearly gain attention in the context of honey bee caste development in general, not

534

only in ovary phenotype determination.

TE D

M AN U

SC

RI PT

524

With respect to the worker ovary markers, oclp-1 and 15-pgdh have as a

536

common feature an expression peak in L5F1 ovaries, and this is also shared with the

537

sdr gene. In contrast, apoLp-III showed a more continuous increase in transcript

538

levels, and thus seems to be of importance or related with later processes of worker

539

ovary development. While generally being highly expressed in the fat body and

540

having a protein function associated with lipid transport in the hemolymph, its

541

differential expression in the developing ovaries is possibly associated with other

542

functions. For instance, in adult honey bees, apoLp-III has been shown to be down-

543

regulated in response to an immune system challenge (Lourenço et al., 2009).

AC C

EP

535

22

ACCEPTED MANUSCRIPT The oclp-1 gene product, which encodes an inhibitor cysteine knot peptide,

545

has been found ubiquitously expressed in all honey bee life cycle stages and tissues,

546

functioning as a conserved antimicrobial or phenoloxidase inhibitor (Bloch and

547

Cohen, 2014). As shown in Tables 2 and 3, several apidaecin and hymenoptaecin

548

peptide genes related to microbial defense functions were indicated as differentially

549

expressed in ovaries, even though this was not validated by RT-qPCR. Such peptides,

550

which have been identified as biologically active in the hemolymph of adult bees,

551

turned out to be inactive in larvae (Casteels et al., 1989). Also, their frequent

552

appearance in differential gene expression analyses of adult and larval honey bees

553

(Flenniken and Andino, 2013; Mao et al., 2015) is possibly related to functions other

554

than microbial defense.

M AN U

SC

RI PT

544

Similarly, the finding that a prostaglandin dehydrogenase could be a consistent

556

marker for worker ovary development may appear quite strange. Nonetheless, even

557

though less studied than in mammals, prostaglandins and eicosanoids are known to

558

play a role in the insect immune response and reproductive biology (Stanley, 2006).

559

In the context of worker ovary development, however, we believe that the 15-PGDH-

560

like protein may act in concert with the SDR protein. Both are short chain

561

dehydrogenase reductases, with 15-PGDH (GB 463670) being a member of the SDR

562

family 36C and SDR (GB54419) belonging to SDR family 11. Insect SDRs are small

563

230-260 aa residue proteins belonging to the NADB-Rossmann superfamily and are

564

characterized by an ADH-short domain. Strikingly, SDRs have appeared in previous

565

studies on differential gene expression related to caste characters in honey bees. In the

566

first one, aiming to identify DEGs in the late larval ovary (Hepperle and Hartfelder,

567

2001), sdr expression was associated with differences in the timing of the ecdysteroid

568

peak in the last instar, and subsequently turned out to be down-regulated by the

AC C

EP

TE D

555

23

ACCEPTED MANUSCRIPT principal honey bee ecdysteroid hormone, makisterone A (Guidugli et al., 2004).

570

Subsequently, and consistent with the current microarray results, this same sdr gene

571

had also been identified as differentially expressed in our previous analysis on larval

572

honey bee ovaries (Humann and Hartfelder, 2011), and more recently we could show

573

that sdr transcripts are overrepresented in optic lobe/retina complexes of worker

574

larvae compared to drones (Marco Antonio and Hartfelder, 2016). The information

575

that now emerged from the current study is that the sdr gene is not only down-

576

regulated by ecdysteroids, but also up-regulated by JH, which makes it an interesting

577

response gene to the caste-specific larval hormone titers (Rachinsky et al., 1990).

SC

RI PT

569

A gene that had caught our attention when seen as overrepresented in worker

579

larval ovaries was hex70b. Hexamerins are conventionally classified as a major class

580

of larval storage proteins (Pan and Telfer, 1996) that are evolutionarily related to

581

prophenoloxidases and arthropod hemocyanins (Burmester, 2002). A general pattern

582

across insect orders appears to be a suppressive effect of JH on hexamerin expression

583

in the late larval fat body (Goodman and Cusson, 2012), and in post-eclosion

584

mosquito females, this repressor effect has been shown to be mediated by the JH

585

receptor methoprene-tolerant (Zou et al., 2013). Furthermore, certain hexamerins have

586

been shown to function as JH binding proteins (Braun and Wyatt, 1996). However, a

587

function that goes beyond transport-related JH binding and which actually sequesters

588

JH and, thus, blocks its action, has so far only been documented in termite caste

589

development (Zhou et al., 2007)

AC C

EP

TE D

M AN U

578

590

In the honey bee, the best studied hexamerin protein is hexamerin 70b, which

591

is not repressed by JH (Cunha et al., 2005). Furthermore, honey bee hexamerins are

592

not only expressed in larval gonads of both sexes (Martins et al., 2010; Martins et al.,

593

2011), but also showed a cytoplasmatic and nuclear localization (Martins and Bitondi,

24

ACCEPTED MANUSCRIPT 2012), colocalizing with a nucleolar protein (Martins and Bitondi, 2016), thus

595

suggesting a regulatory function well beyond what one would expect from a typical

596

storage protein. In a transcriptional profiling study on honey bee queen and worker

597

larvae, hexamerins were also seen as differentially expressed at different time points

598

(Cameron et al., 2013), and in a functional assay these authors found that the RNAi-

599

mediated knockdown of the hex70b gene promoted the development of bees with a

600

queen-like phenotype. Nonetheless, in contrast to their comprehensive transcriptome

601

analysis done on whole body RNA extracts, we found that the hex70b gene is more

602

highly expressed in fourth and early fifth-instar queen ovaries (Figure 2).

603

Furthermore, its expression is down-regulated by JH (Figure 3). These findings

604

indicate that there are tissue-specific responses to JH in the expression of certain

605

hexamerin genes that need to be taken into account in future studies on honey bee

606

caste development.

M AN U

SC

RI PT

594

TE D

607 608

4.4. Defining ovariole numbers – a natural experiment on phenotypic plasticity in

609

honey bee females with implications for drones Together with previous findings, especially of two long noncoding RNAs that

611

are differentially expressed in the larval honey bee ovaries (Humann et al., 2013), the

612

current data provide insights into the molecular architecture underlying the

613

development of the divergent ovary phenotypes in the two female castes. This may

614

lead up to a much wider question, the genetic architecture of the exaggerated gonad

615

morphology in the genus Apis. As stated above, ovariole numbers in this genus by far

616

exceed the basal state in Apoidea (3 ovarioles per ovary) and also in the other

617

members of the Corbiculata clade (four ovarioles per ovary). Strikingly, the same

618

exaggerated morphology is also present in the male sex, the drones, where over 200

AC C

EP

610

25

ACCEPTED MANUSCRIPT testiolar tubules are present in each of the early larval testes (Cruz-Landim, 2009;

620

Zander, 1916). With the exception of honey bee drones, the males of all other bees

621

have only three (basal state) to four (Corbiculata) seminiferous tubules per testis, so

622

the question arises as to whether this trait shared by both sexes in the genus Apis may

623

have a common underlying molecular architecture. In evo-devo terms, gonad

624

development in this genus could thus be interpreted in the framework of the cross-

625

sexual transfer hypothesis (West-Eberhard, 2003), with the additional twist of

626

alternative ovary phenotypes in the female sex generated through massive

627

programmed cell death in the worker larval ovary. As recently stated by Ronai et al.

628

(2015), programmed cell death in the reproductive system is a common theme

629

underlying worker sterility in social Hymenoptera, wherein the massive cell death

630

seen in the larval honey bee ovary is just an extreme and probably recent and lineage-

631

specific version. Their hypothesis is that the incorporation of programmed cell death

632

into a social life style that hinges on worker sterility may have been derived from a

633

general mechanism whereby female insects temporarily block the oogenesis process

634

at specific control points. Obviously, with the current data we are just barely

635

scratching the surface of this much wider evo-devo perspective on reproductive traits

636

in social insects.

SC

M AN U

TE D

EP

AC C

637

RI PT

619

638

Acknowledgements

639

We thank Luis Roberto Aguiar for his help in rearing the many queens needed for this

640

study, Gustavo Jacomini Tiberio for help with the ovary dissections, and Dr. Karina

641

R. Guidugli-Lazzarini for assistance in the RT-qPCR assays. The study received

642

financial support from Fundação de Amparo à Pesquisa do Estado de São Paulo

26

ACCEPTED MANUSCRIPT 643

(FAPESP, grant numbers 2014/08147-3, 2011/03171-5, 2012/01808/-9) and Conselho

644

Nacional de Desenvolvimento Científico e Tecnológico (CNPq, grant number

645

303401/2014-1).

646

RI PT

647

References

649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683

Barchuk, A.R., Cristino, A.S., Kucharski, R., Costa, L.F., Simões, Z.L.P., Maleszka, R., 2007. Molecular determinants of caste differentiation in the highly eusocial honeybee Apis mellifera. BMC Dev Biol 7, e70. Belles, X., Santos, C.G., 2014. The MEKRE93 (Methoprene tolerant-Kruppel homolog 1-E93) pathway in the regulation of insect metamorphosis, and the homology of the pupal stage. Insect Biochem Mol Biol 52, 60-68. Bloch, G., Cohen, M., 2014. The expression and phylogenetics of the Inhibitor Cysteine Knot peptide OCLP1 in the honey bee Apis mellifera. J of Insect Physiol 65, 1-8. Braun, R.P., Wyatt, G.R., 1996. Sequence of the hexameric juvenile hormone-binding protein from the hemolymph of Locusta migratoria. J Biol Chem 271, 3175631762. Büning, J., 1994. The Insect Ovary. Chapman & Hall, London. Burmester, T., 2002. Origin and evolution of arthropod hemocyanins and related proteins. J Comp Physiol B Biochem Mol Biol 172, 95–107. Cai, M.J., Li, X.R., Pei, X.Y., Liu, W., Wang, J.X., Zhao, X.F., 2014. Heat shock protein 90 maintains the stability and function of transcription factor Broad Z7 by interacting with its Broad-Complex-Tramtrack-Bric-a-brac domain. Insect Mol Biol 23, 720-732. Cameron, R.C., Duncan, E.J., Dearden, P.K., 2013. Biased gene expression in early honeybee larval development. BMC Genomics 14, e903. Cardinal, S., Straka, J., Danforth, B.N., 2010. Comprehensive phylogeny of apid bees reveals the evolutionary origins and antiquity of cleptoparasitism. Proc Natl Acad Sci U S A 107, 16207-16211. Casteels, P., Ampe, C., Jacobs, F., Vaeck, M., Tempst, P., 1989. Apidaecins antibacterial peptides from honeybees. Embo J 8, 2387-2391. Cruz-Landim, C., 2009. Abelhas: Morfologia e Função de Sistemas. Editora Unesp, São Paulo. Cruz-Landim, C., Reginato, R.D., Imperatriz-Fonseca, V.L., 1998. Variation on ovariole number in Meliponinae (Hymenoptera, Apidae) queen's ovaries, with comments on ovary development and caste differentiation. Pap. Avul. Zool. (São Paulo) 40, 289-296. Cunha, A.D., Nascimento, A.M., Guidugli, K.R., Simões, Z.L.P., Bitondi, M.M.G., 2005. Molecular cloning and expression of a hexamerin cDNA from the honey bee, Apis mellifera. J Insect Physiol 51, 1135-1147.

AC C

EP

TE D

M AN U

SC

648

27

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Dallacqua, R.P., Bitondi, M.M.G., 2014. Dimorphic ovary differentiation in honeybee (Apis mellifera) larvae involves caste-specific expression of homologs of Ark and Buffy cell death genes. PLoS One 9, e98088. Dedej, S., Hartfelder, K., Aumeier, P., Rosenkranz, P., Engels, W., 1998. Caste determination is a sequential process: effect of larval age at grafting on ovariole number, hind leg size and cephalic volatiles in the honey bee (Apis mellifera carnica). J Apicult Res 37, 183-190. Dietz, A., Hermann, H.R., Blum, M.S., 1979. Role of exogenous JH-I, JH-III and Anti-JH (Precocene-II) on queen induction of 4.5-day-old worker honey bee larvae. J Insect Physiol 25, 503-512. Fittkau, E.J., Kilinge, H., 1973. On biomass and trophic structure of the central Amazonian rain forest ecosystem. Biotropica 5, 2-14. Flenniken, M.L., Andino, R., 2013. Non-specific dsRNA-mediated antiviral response in the honey bee. PLoS One 8, e77263. Gentleman, R.C., Carey, V.J., Bates, D.M., Bolstad, B., Dettling, M., Dudoit, S., Ellis, B., Gautier, L., Ge, Y.C., Gentry, J., Hornik, K., Hothorn, T., Huber, W., Iacus, S., Irizarry, R., Leisch, F., Li, C., Maechler, M., Rossini, A.J., Sawitzki, G., Smith, C., Smyth, G., Tierney, L., Yang, J.Y.H., Zhang, J.H., 2004. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol 5, R80. Goodman, W.G., Cusson, M., 2012 The juvenile hormones, in: Gilbert, L.I. (Ed.), Insect Endocrinology. Academic Press, Oxford, pp. 310-365. Guidugli, K.R., Hepperle, C., Hartfelder, K., 2004. A member of the short-chain dehydrogenase/reductase (SDR) superfamily is a target of the ecdysone response in honey bee (Apis mellifera) caste development. Apidologie 35, 3747. Gutierrez-Beltran, E., Denisenko, T.V., Zhivotovsky, B., Bozhkov, P.V., 2016. Tudor staphylococcal nuclease: biochemistry and functions. Cell Death Diff [E-pub, doi:10.1038/cdd.2016.93]. Hartfelder, K., Emlen, D., 2012. Endocrine control of insect polyphenism, in: Gilbert, L.I. (Ed.), Insect Endocrinology. Academic Press, Oxford, pp. 464-522. Hartfelder, K., Guidugli-Lazzarini, K.R., Cervoni, M.S., Santos, D.E., Humann, F.C., 2015. Old threads make new tapestry - rewiring of signalling pathways underlies caste phenotypic plasticity in the honey bee, Apis mellifera L. . Adv Insect Physiol 48, 1-36. Hartfelder, K., Steinbrück, G., 1997. Germ cell cluster formation and cell death are alternatives in caste-specific differentiation of the larval honey bee ovary. Invertebr Reprod Dev 31, 237-250. He, Q.Y., Wen, D., Jia, Q.Q., Cui, C.L., Wang, J., Palli, S.R., Li, S., 2014. Heat Shock Protein 83 (Hsp83) facilitates Methoprene-tolerant (Met) nuclear import to modulate juvenile hormone signaling. J Biol Chem 289, 27874-27885. Hepperle, C., Hartfelder, K., 2001. Differentially expressed regulatory genes in honey bee caste development. Naturwissenschaften 88, 113-116. Hölldobler, B., Wilson, E.O., 1990. The Ants. Springer, Berlin. Humann, F.C., Hartfelder, K., 2011. Representational Difference Analysis (RDA) reveals differential expression of conserved as well as novel genes during castespecific development of the honey bee (Apis mellifera L.) ovary. Insect Biochem Mol Biol 41, 602-612.

AC C

684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731

28

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Humann, F.C., Tiberio, G.J., Hartfelder, K., 2013. Sequence and expression characteristics of long noncoding RNAs in honey bee caste development potential novel regulators for transgressive ovary size. PLoS One 8, e78915. Iwata, K., 1955. The comparative anatomy of ovary in Hymenoptera. Part I. Mushi 29, 17-34. Jindra, M., Palli, S.R., Riddiford, L.M., 2013. The juvenile hormone signaling pathway in insect development. Annu Rev Entomol 58, 181-204. Kucharski, R., Maleszka, J., Foret, S., Maleszka, R., 2008. Nutritional control of reproductive status in honeybees via DNA methylation. Science 319, 18271830. Leimar, O., Hartfelder, K., Laubichler, M.D., Page, R.E., 2012. Development and evolution of caste dimorphism in honeybees - a modeling approach. Ecol Evol 2, 3098-3109. Linksvayer, T.A., Kaftanoglu, O., Akyol, E., Blatch, S., Amdam, G.V., Page, R.E., 2011. Larval and nurse worker control of developmental plasticity and the evolution of honey bee queen-worker dimorphism. J Evol Biol 24, 1939-1948. Linksvayer, T.A., Rueppell, O., Siegel, A., Kaftanoglu, O., Page, R.E., Amdam, G.V., 2009. The genetic basis of transgressive ovary size in honeybee workers. Genetics 183, 693-707. Liu, Z.G., Ji, T., Yin, L., Shen, J., Shen, F., Chen, G.H., 2013. Transcriptome sequencing analysis reveals the regulation of the hypopharyngeal glands in the honey bee, Apis mellifera carnica Pollmann. PLoS One 8, e81001. Lourenço, A.P., Mackert, A., Cristino, A.D., Simões, Z.L.P., 2008. Validation of reference genes for gene expression studies in the honey bee, Apis mellifera, by quantitative real-time RT-PCR. Apidologie 39, 372-U333. Lourenço, A.P., Martins, J.R., Bitondi, M.M.G., Simões, Z.L.P., 2009. Trade-off between immune stimulation and expression of storage protein genes. Arch Insect Biochem 71, 70-87. Makert, G.R., Paxton, R.J., Hartfelder, K., 2006. Ovariole number - a predictor of differential reproductive success among worker subfamilies in queenless honeybee (Apis mellifera L.) colonies. Behav Ecol Sociobiol 60, 815-825. Mao, W., Schuler, M.A., Berenbaum, M.R., 2015. A dietary phytochemical alters caste-associated gene expression in honey bees. Sci Adv 1, e1500795. Marco Antonio, D.S., Hartfelder, K., 2016. Toward an understanding of divergent compound eye development in drones and workers of the honeybee (Apis mellifera L.): a correlative analysis of morphology and gene expression. J Exp Zool Mol Dev Evol, doi: 10.1002/jez.b.22696. Martins, G.F., Serrão, J.E., 2004. Changes in the reproductive tract of Melipona quadrifasciata anthidioides (Hymenoptera: Apidae, Meliponini) queen after mating. Sociobiology 44, 241-254. Martins, J.R., Anhezini, L., Dallacqua, R.P., Simões, Z.L.P., Bitondi, M.M.G., 2011. A honey bee hexamerin, HEX 70a, is likely to play an intranuclear role in developing and mature ovarioles and testioles. PLoS One 6, e29006. Martins, J.R., Bitondi, M.M.G., 2012. Nuclear immunolocalization of hexamerins in the fat body of metamorphosing honey bees. Insects 3, 1039-1055. Martins, J.R., Bitondi, M.M.G., 2016. The HEX 110 hexamerin is a cytoplasmic and nucleolar protein in the ovaries of Apis mellifera. PLoS One 11, e0151035. Martins, J.R., Nunes, F.M.F., Cristino, A.S., Simões, Z.P., Bitondi, M.M.G., 2010. The four hexamerin genes in the honey bee: structure, molecular evolution and

AC C

732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780

29

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

function deduced from expression patterns in queens, workers and drones. BMC Mol Biol 11, e23. Michelette, E.R.D., Soares, A.E.E., 1993. Characterization of preimaginal developmental stages in Africanized honey bee workers (Apis mellifera L). Apidologie 24, 431-440. Michener, C.D., 2007. The Bees of the World, 2nd ed. The Johns Hopkins University Press, Baltimore, MD. Moda, L.M., Vieira, J., Freire, A.C.G., Bonatti, V., Bomtorin, A.D., Barchuk, A.R., Simões, Z.L.P., 2013. Nutritionally driven differential gene expression leads to heterochronic brain development in honeybee aastes. PLoS One 8, e64815. Niu, D., Zheng, H., Corona, M., Lu, Y., Chen, X., Cao, L., Sohr, A., Hu, F., 2014. Transcriptome comparison between inactivated and activated ovaries of the honey bee Apis mellifera L. Insect Mol Biol 23, 668-681. Pan, M.L., Telfer, W.H., 1996. Methionine-rich hexamerin and arylphorin as precursor reservoirs for reproduction and metamorphosis in female luna moths. Arch Insect Biochem 33, 149-162. Rachinsky, A., Strambi, C., Strambi, A., Hartfelder, K., 1990. Caste and metamorphosis - hemolymph titers of juvenile hormone and ecdysteroids in last instar honeybee larvae. Gen Comp Endocrinol 79, 31-38. Reginato, R.D., Cruz-Landim, C., 2001. Differentiation of the worker's ovary in Apis mellifera L. (Hymenoptera, Apidae) during life of the larvae. Invertebr Reprod Dev 39, 127-134. Rembold, H., 1987. Caste specific modulation of juvenile hormone titers in Apis mellifera. Insect Biochem 17, 1003-1006. Rembold, H., Czoppelt, C., Rao, P.J., 1974. Effect of juvenile hormone treatment on caste differentiation in honeybee, Apis mellifera. J Insect Physiol 20, 11931202. Ronai, I., Vergoz, V., Oldroyd, B.P., 2015. The mechanistic, genetic, and evolutionary basis of worker sterility in the social Hymenoptera. Adv Insect Physiol 48, 251-317. Rozen, J.G., 1986. Survey of the number of ovarioles in various taxa of bees (Hymenoptera, Apoidea). Proc Entomol Soc Wash 88, 707-710. Schmidt Capella, I.C., Hartfelder, K., 1998. Juvenile hormone effect on DNA synthesis and apoptosis in caste-specific differentiation of the larval honey bee (Apis mellifera L.) ovary. J Insect Physiol 44, 385-391. Snodgrass, R.E., 1956. Anatomy of the Honey Bee. Cornell University Press, Ithaca. Stanley, D., 2006. Prostaglandins and other eicosanoids in insects: biological significance. Annu Rev Entomol 51, 25-44. Thompson, G.J., Kucharski, R., Maleszka, R., Oldroyd, B.P., 2006. Towards a molecular definition of worker sterility: differential gene expression and reproductive plasticity in honey bees. Insect Mol Biol 15, 637-644. West-Eberhard, M.J., 2003. Developmental Plasticity and Evolution. Oxford University Press, Oxford. Wheeler, W.M., 1910. The Ants: their structure, develoment and behavior. Columbia University Press, New York. Whitfield, C.W., Cziko, A.M., Robinson, G.E., 2003. Gene expression profiles in the brain predict behavior in individual honey bees. Science 302, 296-299. Zander, E., 1916. Die Ausbildung des Geschlechts bei der Honigbiene (Apis mellifica L.). Z angew Entomol 3, 1-74.

AC C

781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829

30

ACCEPTED MANUSCRIPT 830 831 832 833 834 835 836

Zhou, X.G., Tarver, M.R., Scharf, M.E., 2007. Hexamerin-based regulation of juvenile hormone-dependent gene expression underlies phenotypic plasticity in a social insect. Development 134, 601-610. Zou, Z., Saha, T.T., Roy, S., Shin, S.W., Backman, T.W.H., Girke, T., White, K.P., Raikhel, A.P., (2013) Juvenile hormone and its receptor, methoprene-torelant, control the dynamics of mosquito gene expression. Proc Natl Acad Sci U S A 110, E2173-E2181

RI PT

837

838

SC

839

AC C

EP

TE D

M AN U

840

31

ACCEPTED MANUSCRIPT 841

FIGURE LEGENDS

842 Figure 1: Microarray results on differentially expressed genes (DEGs) in honey bee

844

ovary development. A) Heatmap representing expression levels of 30 DEGs

845

(Benjamini-Hochberg adjusted P-value < 0.05) in fourth (L4) and early fifth-instar

846

(L5) queen (Q) and worker (W) larval ovaries; the extensions r and g in the sample

847

names of the X-axis refer to red and green channel, respectively. These were entered

848

separately in the data analysis to increase statistical power. The gene denominators on

849

the Y-axis are the spot IDs of the microarray. B) Biological function categories of

850

DEGs with known or predicted function. For gene names and lists see Tables 2 and 3.

M AN U

SC

RI PT

843

851

Figure 2: Temporal expression profiles of 15 genes selected for real-time quantitative

853

PCR analysis. The genes were chosen based on differential expression level in the

854

microarray analysis (log2FC in Tables 2 and 3) and biological function. They are

855

organized as genes coding for enzymes (A), transcription factors or signaling pathway

856

proteins (B), storage proteins (C), immune defense factors (D), or other (unknown)

857

functions (E). The RT-qPCR analyses were done on RNA extracted from pools of 10

858

pairs of queen and worker ovaries each dissected from fourth instar (L4), early

859

(L5F1), mid (L5F2) and late feeding stage (L5F3) fifth-instar larvae. Relative

860

expression was calculated as 2-∆∆Ct, by normalization against the endogenous control

861

genes rp49 and actin and calibration against one of the worker L4 ovary samples.

862

Data points represent means ± S.E.M, for three independent biological samples per

863

caste and developmental stage. Statistical significance was analyzed by two-way

864

ANOVA followed by Holm-Sidak multiple comparison tests; * p<0.05.

AC C

EP

TE D

852

865

32

ACCEPTED MANUSCRIPT Figure 3: Effect of juvenile hormone (JH) treatment genes confirmed by prior RT-

867

qPCR time course analysis as differentially expressed in queen and worker ovaries.

868

L4 worker larvae received a topical application of 10 µg JH III in acetone and were

869

collected 6 h later for ovary dissection. Control larvae received an equivalent amount

870

of solvent or were not treated. RT-qPCR analyses were done on three independent

871

biological samples of ovary pools. Bars represent means ± S.E.M. Statistical

872

significance was analyzed by Kruskal-Wallis tests; * p<0.05.

AC C

EP

TE D

M AN U

SC

873

RI PT

866

33

ACCEPTED MANUSCRIPT 874

Supplementary material

875 Figure S1: Correlation plot on log2FC values for the genes selected from the

877

microarray analysis that were further analyzed by RT-qPCR. The RT-qPCR log2FC

878

values shown are the means calculated from the three substages of the fifth larval

879

instar. Linear regression: r2=0.3492, F1,12=6.483, p=0.261; Spearman rank correlation:

880

r=0.6088, two-tailed p=0.0209.

RI PT

876

SC

881

Table S1: Two-way ANOVA results for the 15 genes selected for RT-qPCR analysis.

883

Statistically significant terms are highlighted on gray background.

AC C

EP

TE D

M AN U

882

34

ACCEPTED MANUSCRIPT

BeeBase ID

GenBank ID

Name

Forward primer (5’-3’)

Reverse primer (5’-3’)

Q Q W W W W I I I I I I

GB46223 GB54419 GB43690 GB54372 GB47546 GB46366 GB55452 GB55030 GB48109 GB50902 GB51697 GB42720 GB41845 GB40976 GB50048

NM_001040223.1 NM_001011620.1 XM_001123053.2 XM_392899.5 NM_001011613.1 XM_623815.2 NM_001114198.1 XM_001120252.2 XM_392686.3 XM_393605.5 NM_001011600.1 XM_624891.4 XM_006564514.1 NM_001160064.1 NM_001011629.1

obp14 sdr hsp60 apid14 15-pgdh apoLp-III oclp-1 acpep1 gpdh hex70b hel25e mapk-3 hsp90 e93

CTGGCATTGATCAGCAAAAAGC CTGGAGCTAACTCGGGCATT GACTACGTGCCGCCTAAAGT GGAGGCGGTACTGCTCTTTT CGGCACGAGAGAATTGGTGT TGTGCCCTGGTGTTACAACT CTCCCGAATTGGAAAGATCA CGAACATCGTTTCAGCTGCC TCGCCATGTACTGGGTGAAC CTGCACAGACCCGAGTGAAT GGGTGACACAGCTGACATGA TTGTTGGTACTCCAGGTCGC TCCCGAGTGGTGTGAAGGTA TGGCAAACAGTTGGTCTCTG GGTGGACGCGTGGATTTGA

GACTGCTTTGATTCCTTGTGGT CGTTTTGGTTGGAAAGATCGCA AAGAGCAGCCTTCACAGCAT TAGCATCCACGCCTGCATTT AGTAGGCGGATCTAGGTTGGT GCGTTTGCGGGATGTCTTTT GCTCAGAGATTTCGCAGCTT AGCTCCTCCTCTGTGATCCT TCCTTTGGTACCTGTGCGAC CAACAACCTGAGCACCGAAC GAGGCCAACATCTTCGGTGA ACGTCCCTACGCATGTCTAA CCTCAGGACCAACGGAATCA ACAGCATGGAGAATCAACTAACCT CGATCGAGACACCGAGAGGA

M AN U

TE D EP

Q W W

SC

L5

AC C

L4

RI PT

Table 1: Genes selected from the microarray gene list for RT-qPCR analysis, their library origin (fourth or fifth larval instar), gene number IDs, provisional gene name abbreviation and primer sequences.

ACCEPTED MANUSCRIPT

Table 2: Genes identified as higher expressed in ovaries of queen larvae. Shown are instar, microarray and BeeBase IDs, predicted gene names, log2 fold change and B statistics. Results are in descending order of log2 fold change. Differentially expressed genes in fourth instar ovaries are on gray background. Entries in bold represent genes that had their transcript levels further analyzed by quantitative PCR. Beebase ID

log2FC

B

AM03163 AM06387

GB46223 -

odorant binding protein 14 (obp14) SET: OGS; Dmel Ortholog - : FBgn0013678

1,045 1,038

0,55 0,60

L5

AM11525

GB54372

heat shock protein 60 (hsp60)

2,541

5,27

L5 L5

AM05659 AM12887R

GB54419 -

short chain dehydrogenase reductase (sdr) [Arrayset: JDE_EST jdeH07_2DEF2, Group10.13_52060_52122]

2,003 1,899

2,62 4,23

L5 L5

AM00359 AM01778

GB47546 BB170022A10A06.5

apidaecin 14 (apid14) EST Bee Brain Normalized/Subtracted Library, BB17

1,542 1,539

3,12 2,63

L5

AM01148

BB170024B20F01.5

EST Bee Brain Normalized/Subtracted Library, BB17

1,523

2,38

L4 L4

predicted gene name

AM02500

GB50048

transcription factor mblk-1 (Mblk-1)

AM07353

GB40976

heat shock protein 90-alpha (hsp90)

L5 L5

AM07393 AM11827

GB50902 GB43180

L5 L5

AM03488 AM01045

GB41845 BB170011A20H02.5

L5

AM05474

GB40977

L5

AM12865R

-

L5

AM00352

GB47546

L5 L5

AM05822 AM00357

GB42720 GB47546

L5

AM00356

GB51306

1,63 0,95

glyceraldehyde-3-phosphate dehydrogenase (gpdh) minor histocompatibility antigen H13-like

1,400 1,385

0,98 2,51

mitogen-activated protein kinase phosphatase-3 (mapk-3) EST Bee Brain Normalized/Subtracted Library, BB17

1,330 1,328

0,80 1,64

tudor-SN

1,316

0,98

1,281

0,92

apidaecin 14 (apid14)

1,257

2,00

helicase at 25E (hel25e) apidaecin 14 (apid14)

1,242 1,229

1,33 1,65

apidaecin (apid73) complement component 1 Q subcomponent-binding protein, mitochondrial-like

1,228

1,94

EP

TE D

[Arrayset: JDE_EST jdeH07_2DEF2, Group10.13_52060_52122]

L5

AM02735

1,179

1,28

L5

AM00354

NM_001011642.1

apidaecin type 22 (apid22)

1,162

1,58

L5

AM00351

NM_001011642.1

apidaecin type 22 (apid22)

1,115

0,87

L5

AM07489

SET: OGS; Dmel Ortholog - C. Elsik Feb_2007: FBgn0038964 akirin

1,094

1,14

1,092

1,05

head Apis mellifera cDNA clone BH10068H09 3

1,086

1,16

deoxyuridine 5-triphosphate nucleotidohydrolase (dUTPase) actin related protein 1 (arp1)

1,083

0,90

1,050

0,86

AM02286 Unknown [Arrayset: UI_EST BI514555]

1,040

0,70

apidaecin type 22 (apid22)

1,039

0,67

hymenoptaecin

1,029

0,83

AM12836R Unknown [Arrayset: JDE_EST jdeC15_]

1,016

0,70

apidaecin 14 (apid14)

1,002

0,54

-

L5

AM06889

GB51067

L5

AM00831

DB774043

L5

AM01471

-

L5

pos2

GB44311

L5

AM02286

L5

AM00355

L5

AM10109

L5

AM12836

L5

AM00358

NM_001011642.1 GB51223 GB47546

AC C

GB44693

1,456 1,414

M AN U

L5 L5

RI PT

microarray ID

SC

instar

ACCEPTED MANUSCRIPT

Table 3: Genes identified as higher expressed in ovaries of worker larvae. Shown are instar, microarray and BeeBase IDs, predicted gene names, log2 fold change and B statistics. Results are in descending order of log2 fold change. Differentially expressed genes in fourth instar ovaries are on gray background. Entries in bold represent genes that had their transcript levels further analyzed by quantitative PCR. instar

microarray ID

Beebase ID

predicted gene name

log2FC

B

AM07204

GB54419

short-chain dehydrogenase/reductase (sdr)

4,084

2,969

L4 L4

AM06564 AM06212

GB43690 GB51696

CG4409-PA (LOC727344) hexamerin 70c (hex70c)

2,518 2,449

3,027 3,102

L4

AM01327

BB170032B10H06.5

EST Bee Brain Normalized/Subtracted Library, BB17

1,892

1,746

L4

AM01720

BI946450.1

bEST46.3' Honeybee brain cDNA library

1,630

2,189

L4

AM12859

GB17785

1,414

0,587

L4

AM12845

BB160015B10B06.5

1,037

0,636

L5

AM07204

GB40137

5,136

5,996

L5 L5

AM11298 AM00326

GB46367 GB30416

L5

AM12664 AM11850

SC

EST Bee Brain Normalized Library, BB16 serine protease 17 (sp17)

3,254 2,260

5,487 5,147

GB55452

15-hydroxyprostaglandin dehydrogenase (15-pgdh) cuticular protein 65Av apolipophorin-III-like protein (apoLp-III)

1,801

4,162

GB55030

hypothetical inhibitor function cysteine knot peptide (oclp-1)

1,554

2,936

1,543 1,505

3,043 3,019

M AN U

L5

Nedd8-like

RI PT

L4

AM03896 AM07163

GB48109 GB45943

L5

AM04995

GB47542

Adam eukaryotic translation initiation factor 3 subunit J

1,505

2,452

L5 L5

AM03494 AM06212

GB51697 GB51696

hexamerin 70b (hex70b) hexamerin 70c (hex70c)

1,444 1,441

2,851 2,442

L5

AM03359

GB44777

CG17127-PA LOC725882

1,433

2,833

L5

AM02239

BB160014B10G10.5

EST Bee Brain Normalized Library, BB16

1,427

2,807

L5 L5

AM03163 AM05331

GB46223 GB43689

odorant binding protein 14 (obp14) CG4409-PA LOC727522

1,418 1,403

1,200 2,381

L5

AM07200

GB40137

1,401

2,350

L5

AM11377

GB52324

serine protease 17 (sp17) chemosensory protein 3 (csp3)

1,383

2,531

L5 L5

AM06564 AM07609

GB43690 GB50121

CG4409-PA (LOC727344) Amci chymotrypsin inhibitor

1,371 1,333

0,744 1,500

L5

AM01191

BB170008A10H01.5

EST Bee Brain Normalized/Subtracted Library, BB17

1,285

1,237

L5

AM02526

BB170027B10A03.5

EST Bee Brain Normalized/Subtracted Library, BB17

1,235

0,505

L5

AM02017

BB160008B20E05.5

EST Bee Brain Normalized Library, BB16

1,140

1,082

AC C

EP

TE D

L5 L5

retinoid-inducible serine carboxypeptidase-like (acpep1) collagen alpha-1(IV) chain-like

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT Highlights

• • •

AC C

EP

TE D

M AN U

SC

•

Molecular signatures of caste-specific ovary development in honey bee larvae were investigated A microarray analysis showed that differential gene expression increases at the beginning of the 5th larval instar Expression differences seen in microarrays were confirmed by RT-qPCR for nine genes, seven also responding to JH An hsp90 gene emerged as a consistent marker for queen ovary development A gene encoding a short chain dehydrogenase reductase is strongly related to the worker ovary phenotype

RI PT

•