Integration of lncRNA–miRNA–mRNA reveals novel insights into oviposition regulation in honey bees

Background The honey bee (Apis mellifera) is a highly diverse species commonly used for honey production and pollination services. The oviposition of the honey bee queen affects the development and overall performance of the colony. To investigate the ovary activation and oviposition processes on a molecular level, a genome-wide analysis of lncRNAs, miRNAs and mRNA expression in the ovaries of the queens was performed to screen for differentially expressed coding and noncoding RNAs. Further analysis identified relevant candidate genes or RNAs. Results The analysis of the RNA profiles in different oviposition phase of the queens revealed that 740 lncRNAs, 81 miRNAs and 5,481 mRNAs were differently expressed during the ovary activation; 88 lncRNAs, 13 miRNAs and 338 mRNAs were differently expressed during the oviposition inhibition process; and finally, 100 lncRNAs, four miRNAs and 497 mRNAs were differently expressed during the oviposition recovery process. In addition, functional annotation of differentially expressed RNAs revealed several pathways that are closely related to oviposition, including hippo, MAPK, notch, Wnt, mTOR, TGF-beta and FoxO signaling pathways. Furthermore, in the QTL region for ovary size, 73 differentially expressed genes and 14 differentially expressed lncRNAs were located, which are considered as candidate genes affecting ovary size and oviposition. Moreover, a core set of genes served as bridges among different miRNAs were identified through the integrated analysis of lncRNA-miRNA-mRNA network. Conclusion The observed dramatic expression changes of coding and noncoding RNAs suggest that they may play a critical role in honey bee queens’ oviposition. The identified candidate genes for oviposition activation and regulation could serve as a resource for further studies of genetic markers of oviposition in honey bees.

152 Furthermore, the expression of each RNA type was analyzed with unsupervised hierarchical 153 clustering with the R package of "pheatmap". To do unsupervised hierarchical clustering, firstly, 154 the expression of RNA was normalized. For normalization of lncRNA and mRNA, the following 155 formula was used: . For normalization of miRNA, the following FPKm = log 10 FPKM + 1 156 formula was applied: (TPM, transcripts per kilobase million). Then the TPm = log 10 TPM + 1 157 euclidean distance was used to measure the degree of similarity between the expression profiles 158 of samples. The method in the package to cluster distance is "complete".
159 Prediction of lncRNA and miRNA target genes 160 The potential trans role of lncRNAs (acting on non-neighboring genes) can be assessed by 161 correlating expression levels between lncRNAs and mRNAs. The trans role of lncRNAs in 162 coding genes was examined based on the expression correlation coefficient (Pearson correlation 163 ≥ 0.95 or ≤ -0.95). To predict miRNAs targets, we searched for the targets in the 3'UTR of genes 164 models. For genes lacking a predicted 3'UTR, the region 1000bp downstream of the stop codon 165 were included. The prediction was performed by Miranda with the following parameter: free 166 energy < -10 kcal/mol and score > 140 (Enright et al. 2003).

GO and Pathway enrichment analysis
228 Functional annotation analysis of target genes of the DE lncRNAs, miRNA and mRNA was 229 performed to identify GO terms and KEGG pathways with higher confidence (Supplemental 230 Table S3, S4 and S5). Because GO terms and pathways enriched with the DE lncRNAs, miRNA 231 and mRNAs were similar to each other, here we only describe the enrichment results of DE 232 mRNAs. In the ovary activation process, most of the enriched GO_BP terms of DE mRNAs 233 were involved in tissue development, energy producing and hormone biosynthesis and 234 metabolism, such as oocyte microtubule cytoskeleton polarization, fatty acid oxidation, 235 neurotrophin signaling pathway, ecdysteroid catabolic process (Supplemental Table S3). In the 236 oviposition inhibition process, contrary to the ovary activation process, several GO terms were 237 not enriched, but enrichment occurred again when oviposition recovered, such as cellular 238 response to transforming growth factor beta stimulus, positive regulation of cyclase activity, 239 post-embryonic hemopoiesis, larval lymph gland hemopoiesis, eye pigment biosynthetic process, 240 and compound eye cone cell fate commitment (Supplemental Table S3). 241 DE mRNAs enrichment (p < 0.05) was seen in KEGG pathways (Supplemental Table S3). Manuscript to be reviewed 246 pathway (Table 2). 247 Chromosomal localization of DE lncRNAs and mRNAs in QTL region for ovary size 248 If the differentially expressed lncRNAs and mRNAs were found located within the confidence 249 interval of the QTL for ovary size, they could be regarded as candidate genes for ovary size and 250 potential candidate genes for oviposition. In this way, 73 candidate genes and 14 lncRNAs 251 (Supplemental Table S6) were identified.
252 Construction of the lncRNA-miRNA-mRNA network 253 The bioinformatic analysis predicted that 469 lncRNAs were targeted by 69 miRNAs and 117 254 lncRNAs acted as decoys to 31 miRNAs. The transcriptome network was constructed based on 255 the lncRNA-miRNA and the miRNA-mRNA relationship pairs. The resulting network consists 256 of 229 lncRNA-miRNA pairs and 225 miRNA-mRNA pairs (Supplemental Fig.S1 and Table S7).
257 To further investigate the potential candidate genes and RNAs for ovary activation and 258 oviposition, a reproductive associated network was constructed containing the DE miRNAs and 259 mRNAs which played specific or suspected roles in reproduction, and the DE lncRNAs and  Table S8).
263 Validation of RNA-Seq data by real time PCR 264 In order to validate the sequencing results, the expression of 5 lncRNAs, 5 mRNAs and 5 265 miRNAs were tested by using real time PCR with the same RNA samples used for sequencing Table S1). The expression profiles of these genes/RNAs detected by real time 267 PCR were consistent with those obtained by sequencing, which confirmed the reliability of our 268 sequencing results. Manuscript to be reviewed 270 In the present study, dynamical lncRNAs, mRNAs and miRNAs expression profiles in ovary 271 activation and oviposition processes in honey bees were identified. However, the complex 272 molecular mechanism behind the oviposition activation and regulation still needs to be illustrated.
273 Representative enriched pathways 274 The gene function analysis showed that DE RNAs enrichment was seen in a number of pathways 275 in ovary activation and/or oviposition regulation process. Some of the pathways are particularly 276 interesting, such as Wnt, hippo, TGF-beta, notch, MAPK, FoxO and mTOR signaling pathways 277 (Fig.3). More than 50% of the genes in those pathways were differently expressed according to 278 283 The Wnt signaling pathway was found to be involved in the development of reproductive system 284 such as the development of ovarian follicules, ovulation and luteinization (Sun & Wang 2003). 285 The hippo signaling pathway was also reported to be related to the regulation of mouse ovarian 286 functional remodeling (Ye et al. 2017). Moreover, the hippo signaling pathway can coordinate 287 with Wnt, TGF-beta and notch signaling pathways affecting organ size in Drosophila (Barry & 288 Camargo 2013). Because after queen mating, the size of ovary will become bigger than the 289 virgin's (Rinderer 1987), we also observed many genes in Wnt, TGF-beta, hippo and notch 290 signaling pathways that were differentially expressed in mated queens compared with virgin 291 queens. It indicated that those pathways may participate in ovarian function remodeling after 292 mating to prepare for oviposition in honey bees. The oocyte growth and development is crucial 308 Further, the studied pathway maps were looked up in KEGG database to assess whether there is 309 a relationship among them. The results showed that they were closely interacting with each other 310 as shown in Fig.3, whereby for example TGF-beta signaling pathway was part of hippo and Wnt 311 signaling pathways. These pathways were enriched both in oviposition activation and oviposition 312 process. Considering roles of these pathways in ovarian function remodeling, oocyte growth and 313 development and other related processes, they are critical for a successful oviposition by 314 complex fine-tuning relationships. 315 Among the DE genes in those pathways, several genes were found to participate in more than 316 one pathway. The gene nejire (Nej, also known as CREB-binding protein (CBP)) participated in 317 three pathways, namely notch, FoxO and TGF-beta signaling pathways. Additionally, Nej was 318 significantly up-regulated in the egg-laying queens compared to virgin queens. Also studies in 319 Drosophlia melanogaster found that Nej was involved in regulation of many pathways during 320 embryo development, through hedgehog, wingless and TGF-beta signaling pathways 321 (Fernandez-Nicolas & Belles 2016). Taken together, we could conclude that Nej may participate 322 in embryonic development in honey bees through notch, FoxO and TGF-beta signaling pathways, 323 and can be considered as the potential candidate genes for oviposition. 324 Genes and lncRNAs co-localized in QTL region for ovary size 325 We compared the location of the DE genes and DE lncRNAs on the honey bee genome available 326 at the NCBI database (Amel_4.5) with one QTL for ovary size. 73 genes and 14 lncRNAs were 327 identified, and some of them together with their key function will be explained further. 328 Among the 73 genes, G2/mitotic-specific cyclin-B3 (CycB3) is the one we paid special attention.  Table 3 and Table 4 show that 31 mRNAs and 36 miRNAs were significantly regulated in caste 358 determination or other reproductive related process, which indicated that they have known or 359 suspected roles of in ovary activation. The oviposition status is positively correlated with two 360 genes (heat shock protein 90 (Hsp90) and Usp) and negatively correlated with ten genes. Hsp90 361 has been reported as a candidate marker gene for caste-specific ovary development (Lago et al. Manuscript to be reviewed 426 The authors thank Wei Feng, who is the beekeeper in the apiary, for his assistance in beekeeping. 427 Supplemental list 428 Figure S1 The lncRNA-miRNA-mRNA network with at least one DE RNA in a lncRNA-465 miRNA pair or a miRNA-mRNA pair.   V, group of virgin queens (n=3); Q, group of egg-laying queens (n=3); C, group of egg-laying inhibited queens (n=3); R, group of egg-laying recovery queens (n=3).

Figure 2(on next page)
The reproductive associated lncRNA-miRNA-mRNA network. Manuscript to be reviewed