Virus and dsRNA-triggered transcriptional responses reveal key components of honey bee antiviral defense

Recent high annual losses of honey bee colonies are associated with many factors, including RNA virus infections. Honey bee antiviral responses include RNA interference and immune pathway activation, but their relative roles in antiviral defense are not well understood. To better characterize the mechanism(s) of honey bee antiviral defense, bees were infected with a model virus in the presence or absence of dsRNA, a virus associated molecular pattern. Regardless of sequence specificity, dsRNA reduced virus abundance. We utilized next generation sequencing to examine transcriptional responses triggered by virus and dsRNA at three time-points post-infection. Hundreds of genes exhibited differential expression in response to co-treatment of dsRNA and virus. Virus-infected bees had greater expression of genes involved in RNAi, Toll, Imd, and JAK-STAT pathways, but the majority of differentially expressed genes are not well characterized. To confirm the virus limiting role of two genes, including the well-characterized gene, dicer, and a probable uncharacterized cyclin dependent kinase in honey bees, we utilized RNAi to reduce their expression in vivo and determined that virus abundance increased, supporting their involvement in antiviral defense. Together, these results further our understanding of honey bee antiviral defense, particularly the role of a non-sequence specific dsRNA-mediated antiviral pathway.

ImageJ was utilized to quantify the pixel count in each band, and the GFP: β-actin ratio was calculated for each sample and plotted on a box-and-whisker; p-values were assessed via Wilcoxon rank-sum tests. The median GFP to actin ratio was 1.0 for virus-infected bees, whereas the median ratio in bees that were co-injected with sp-dsRNA was 0.3 (p = 0.03). Similarly, co-injection of virus and ns-dsRNA reduced the relative production of SINV-GFP as compared to co-injected with virus and NTPs, with median ratios of 0. 5   The abdomens of bees (A) 48 and (B) 72 hours post-infection (hpi) have reduced relative virus abundance when treated with dsRNA, whether it is sequence-specific to the virus (sp-dsRNA, dotted purple) or nonspecific (ns-dsRNA, checkered blue). This result was consistent in three different biological replicates, including bees from different colonies. Treatment with NTPs (wavy orange lines) did not result in different relative virus abundance. Biological replicate 1 (rep 1) includes bees that were utilized for RNAseq analysis for which their relative virus abundance is also shown in Fig 2. Biological replicates 2 and 3 were bees collected from different colonies and likely represent bees with different genetic compositions. Relative virus abundance was assessed using qPCR and ΔΔCT analysis, using Am rpl8 as the house keeping gene; statistical differences between treatment and virus-infected bees (n=10) were performed using student's t-test, *p< 0.05, **p< 0.005. The bars represent standard error of the mean. Venn diagrams were utilized to identify shared and unique differentially expressed genes (DEGs) among treatment groups within each time point post-infection. By identifying DEGs that were shared by virus-infected bees (green), dsRNA-treated bees (orange), bees treated with both virus and virus-specific dsRNA (sp-dsRNA, blue), and bees treated with virus and nonspecific dsRNA (ns-dsRNA, purple), we could identify the number of genes that are induced by biologically relevant levels of dsRNA at (A) 6 hpi, (B) 48 hpi, and (C) 72 hpi, as well as identify DEGs specific to each treatment or each specific comparison of treatments. (D) The shared DEGs in all virus-infected bees at 48 hpi and 72 hpi were too large to visualize in Venn diagrams, therefore the number of shared DEGs were expressed in table format. In order to identify genes that may contribute to enhanced antiviral defense in dsRNA-treated bees (Fig 2), the shared genes with increased expression in all virus-infected and dsRNA treated bees, but not shared with the bees that were only infected with virus (five genes), were assessed. Five genes exhibited increased expression in virus-infected and dsRNA treated bees, but not in bees only infected with virus. There were two genes with decreased expression in all virus and dsRNA treated bees, but not in virus treated bees. Lists of the specific shared genes in all Venn diagrams (A-F) are provided in Tables S9-S13  Transcriptome level sequencing (RNAseq) identified hundreds of honey bee (Apis mellifera) genes that exhibited differential expression in the context of virus infection and/or dsRNA-treatment, and qPCR was utilized to validate RNAseq results of 10 genes that exhibited increased expression in bees 48 hpi (q < 0.05 after BH correction, Tables S5 and S6). The nine DEGs exhibited increased expression in virus-infected bees (green stripes) and virus and dsRNA treated bees (sp-and ns dsRNA, dotted purple and checkered blue) via qPCR analysis including hsp90, abaecin,probable cyclin-dependent kinase, ago2, dicer, igfn3-10, orb-2-like, mfs-type transporter, and fam102b; whereas jra did not, which was expected based on RNAseq analysis. Two DEGs, jra, and fam102b, exhibited increased expression in the dsRNA (dotted orange) bees, which was expected. Relative gene transcript abundance was assessed using qPCR and ΔΔCT analysis, using Am rpl8 as the house keeping gene; expression in each treatment group was compared to mock-infected controls. Thirteen genes were also qPCR validated in bees at 72 hpi (Fig 4). Statistical differences in gene expression between mock-infected and treatment bees (n=5) were performed using Welch's two sample t-tests, *p< 0.05, **p< 0.005. The bars represent standard error of the mean. Transcriptome analysis (RNAseq) of virus-infected and/or dsRNA-treated bees identified numerous candidate honey bee antiviral genes that exhibited greater expression as compared to mock-infected controls (q-value < 0.05 after BH correction), including the genes (A) hsp90, (B) cyclin-dependent serine/threonine kinase, (C) ago2, (D) dicer, (E) mfs-transporter, (F) formin-j, (G) tet-2, (H) orb-2, (I) 3a1-like, (J) and igfn 3-10. qPCR analysis was performed for these genes on the bees that were sequenced (rep 1). (A-J) In order to confirm that these genes consistently exhibit increased expression in response to virus infection in honey bees with different genetic backgrounds, qPCR was performed on pooled samples (5 bees per pool) from virus-infected bees that were collected from different honey bee colonies (reps 2 and 3). (A-F) Several of these genes displayed increased expression via qPCR in all three biological replicates (A-D, G, J), but not as robustly, perhaps due to the fact that the samples were pooled as opposed to individually analyzed as in rep 1. (E, H, I) Other genes, displayed increased expression in 2 of the 3 biological replicates. (F) One gene, igfn 3-10, only exhibited increased expression in the sequenced bees (rep 1) and not in the other biological replicates (rep 2 and 3). Normalized gene transcript abundance was assessed using qPCR and ΔΔCT analysis using Am rpl8 as the house keeping gene. The bars are standard error of the mean. RNASeq analysis determined that reads aligning to LOC725387 were more abundant in virus-infected bees. To identify the gene or genes encoded by these differentially expressed reads, the consensus nucleotide sequence was used to query the NCBI Nucleotide collection (nr/nt) and A. mellifera databases using blastn, Sanger sequencing was performed to verify transcript sequence and length, and the results were evaluated using Geneious.
(A) The updated LOC725387 A. mellifera probable cyclin-dependent serine/threonine-protein kinase transcript ( A. mellifera To further investigate the role of dicer and probable cyclin-dependent kinase in honey bee antiviral defense, we utilized RNAi-mediated gene knock down to reduce their expression and qPCR to confirm reduced target gene expression and determine its impact on virus abundance (Fig 7). (A) In bees treated with dsRNA in the absence of virus at 48 and 72 hpi, cyclin dependent kinase expression was decreased in kinase specific dsRNA treated bees (spotted orange) by 46% (*p<0.05) and 50% (**p< 0.005) as compared to the ns-dsRNA controls (red horizontal lines). In virus and dsRNA-treated bees at 48 and 72 hpi, cyclin dependent kinase expression was decreased in kinase specific dsRNA treated bees (black diagonal lines) by 40% (*p<0.05) and 30% (*p< 0.05) as compared to the ns-dsRNA controls (checkered blue). (B) In bees treated with dsRNA in the absence of virus at 48 and 72 hpi, dicer expression was decreased in dicer specific dsRNA treated bees (vertical purple lines) by 56% (*p<0.05) and 31% (*p< 0.05) as compared to the ns-dsRNA controls (red horizontal lines). In virus-infected bees at 48 hpi, dicer expression was decreased in dicer specific dsRNA treated bees (spotted blue) by 32% (*p< 0.05), but it was not reduced at 72 hpi as compared to the ns-dsRNA controls. (C) The virus abundance of kinase and dicer specific dsRNA-treated bees had a trend of increased virus abundance compared to the ns-dsRNA controls, but this was not significant at 48 hpi. Virus abundance in kinase and dicer specific dsRNA-treated bees at 72 hpi was increased by 48% (*p< 0.05) and 44% (*p< 0.05) compared to the ns-dsR-NA control. Relative gene transcript abundance was assessed using qPCR and ΔΔCT analysis, using Am rpl8 as the house keeping gene; expression in each treatment group was compared to mock-infected controls. Percent relative virus abundance for each sample was determined via qPCR and ΔΔCT analysis using Am rpl8 as the house keeping gene. Statistical differences between treatment (n=10) and virus-infected and ns-dsRNA bees were determined using one-sided Student's t-tests, *p< 0.05, **p< 0.005. The bars are standard error of the mean.

Sindbis virus (SINV-GFP) infection trials
Currently, there are no infectious honey bee virus clones, and though studies with semipurified honey bee virus preparations have provided valuable information 3-5 , we utilized a recombinant model virus, Sindbis virus expressing green fluorescent protein (SINV-GFP) 1,6 . There are several advantages to utilizing this virus including the ability to control the dose of virus inoculum, monitor the progression of virus infection using GFP, and the assurance that the honey bees were not previously infected with, nor exposed to, SINV-GFP. In addition, Sindbis virus does not encode a suppressor of RNAi (VSR) 7 . We and others have used SINV-GFP to investigate honey bee 1 , fruit fly 6 , and mosquito 8 antiviral defense mechanisms, thus facilitating comparison of immune responses in both natural mosquito hosts and non-native hosts (i.e., honey bee and fruit fly) that have not co-evolved with this virus. In brief, honey bees were immobilized via incubation at 4ºC for 20 minutes and injected in the thorax with 3,750 plaque forming units (PFUs) of SINV-GFP 1 diluted in 2 µl of 10 mM Tris buffer (pH 7.5) using a Harbo large capacity syringe equipped with disposable needles (Honey Bee Insemination Service; http://www.honeybeeinsemination.com/equipment2.html). The needles were prepared from borosilicate capillary tubes (0.8-1.10 x 100 mm) with a micropipette puller (Narishige Model PC-10, East Meadow, New York, USA). To investigate the role of dsRNA in honey bee antiviral defense, SINV-GFP was inoculated with multiple species and lengths of dsRNA (1 µg each), including virus-specific dsRNA (sp-dsRNA, 928 bp), nonspecific dsRNA corresponding to Drosophila C virus sequence (ns-dsRNA, 1,017 bp), and luciferase sequence (LUC dsRNA, 355 bp) (Supplementary Table S1). Bees were also co-injected with 1 µg high molecular weight polyinosinic-polycytidylic acid (poly(I:C)), InvivoGen, San Diego, California, USA), a synthetic mimic of dsRNA, or 1 µg nucleoside triphosphates (NTP), the positive and negative controls, respectively. After injection, bees typically recovered after 5 minutes at room temperature. Bees rarely died after injection (i.e., < 6%) and if so death was attributed to poor injection technique, as it was not associated with the substance injected (e.g., buffer, dsRNA, virus) and those bees were removed from the study. Mock-infection controls were also performed. Bees were collected at 6, 48, or 72 hours post-infection (hpi). This time course was established to examine both early and late immune responses within a time frame that allowed for virus dissemination and infection, while maintaining optimal conditions for bees housed within the laboratory setting 1 . For each experimental treatment two additional biological repetitions that utilized bees from different colonies were evaluated at 48 and 72 hpi.
dsRNA preparation dsRNA was generated by in vitro transcription with T7 RNA polymerase 6,9 . T7 promoter containing dsDNA PCR-products (1-10 µg) were amplified using the primers listed in Supplementary annealed by heating the reaction to 100ºC for 5 minutes and then slowly cooling to room temperature. dsRNA products were purified by phenol:chloroform extraction followed by ethanol precipitation; dsRNA for injection was suspended in 10 mM Tris pH 7.5. dsRNA quality was assessed by agarose gel electrophoresis and spectrophotometer (NanoDrop 2000c, ThermoScientific, Waltham, MA, USA).The dsRNA quantity based on agarose gel band intensity was assessed using ImageJ 10 .

dsRNA-mediated gene knockdown
The expression of two candidate antiviral genes: dicer (XM_016917734.1) and a novel transcript with 91% sequence identity with the Apis cerana probable cyclin-dependent serine/threonine-protein kinase DDB_G0292550 (XM_017051141.1), was reduced by RNAi-mediated gene knockdown (i.e., bees were injected with 1 µg of gene-specific dsRNA) (Supplementary Table S1). In order to assess the effects of gene knockdown on virus abundance, bees were infected with SINV-GFP (using methods as above) and co-injected with either gene-specific or nonspecific (DCV-specific) dsRNA (control).

Fluorescence Microscopy
The abdomens of individual honey bee bees were imaged using a Zeiss Stereoscope Stemi SV11 Apo, equipped with a Jenoptik Camera ProgResC14 Plus (Optical Systems GmBH P-07739 Jena). Fluorescence images were taken under fluorescent light with a GFP filter using standardized camera and exposure settings (i.e, 20x magnification, gain 60 and 150 ms exposure time). White light images were also taken 20x magnification (50 ms exposure).

Honey bee protein lysate preparation and analysis
Bees were dissected into head, thorax, and abdomen. The thoraxes and abdomens of bees collected at 72 hours post-infection were individually homogenized in 400 µl sterile H 2 O with two sterile glass beads (5 mm) via bead beating for 1.5 minutes. The thorax and abdomens of bees collected at 6 and 48 hours post-infection were individually homogenized in 400 µl sterile H 2 O with one sterile 4.5 mm steel ball bearing using a Tissue Lyzer II (Qiagen, Hilden, Germany) for 1 minute. Lysates were clarified by spinning for 12 minutes at 12,000 x g.
For Western blot analysis, protein-containing honey bee lysates were combined with Laemmli buffer (95ºC for 3 min), electrophoresed on 12% acrylamide gels (Mini-PROTEAN TGX, BioRad, Hercules, CA), transferred to Immobilon-P PVDF membranes (EMD Millipore Corporation, Billerica, MA, USA), blocked with 5% milk in TBS buffer containing 0.1% Tween-20, incubated at 4ºC overnight with primary antibodies: either α-GFP (sc-8334; Santa Cruz Biotech, Santa Cruz, CA, USA) or α-β-actin (#4967L; Cell Signaling Technology, Danvers, MA, USA), washed, and then incubated with Horseradish peroxidase-conjugated anti-rabbit secondary antibody (ECL, GE Healthcare) for one hour at room temperature. The Western blots were developed with SuperSignal West Femto Maximum Sensitivity Substrate (Pierce, ThermoScientific, Waltham, MA, USA) and imaged with the Syngene G:Box F3 (software: G:Box Chemi-XR5). ImageJ 10 was used to quantify the sum of the pixels within β-actin and GFP regions in order to calculate and compare SINV-GFP levels in samples. The number of pixels within GFP regions were normalized based on the number of β-actin pixels as a loading control. Treatment groups in each Western blot set were assessed for statistical differences using Wilcoxon ranked-sum test.

Honey bee RNA isolation and purification
TRizol reagent (Invitrogen), 400 µl was added to 400 µl of individual bee thorax and abdomen homogenate, and RNA was isolated according to manufacturer's instructions. Prior to gene expression analysis by RNASeq or qPCR, RNA was further purified using Qiagen RNAeasy columns including on column DNase Treatment (Qiagen) to remove DNA from samples. RNA was quantified using a Thermo Scientific NanoDrop 2000 Spectrophotometer (Waltham, MA, USA).

Quantitative PCR (qPCR)
Quantitative PCR was utilized to examine the relative abundance of virus and honey bee host gene expression in each sample using previously described methods that are in accordance with published guidelines 11 . All qPCRs reactions were performed in triplicate using 2 µL of cDNA as template. Each 20 µl reaction was composed of cDNA template, 1X SYBR Green (Invitrogen, Cat.57563), 1X Choice Taq Master Mix (Denville Scientific Inc., Holliston, MA), 3 mM MgCl 2 , and forward and reverse primers (600 nM each). A CFX Connect Real Time instrument (BioRad, Hercules, CA) was utilized for qPCR, the thermo-profile for virus (e.g., SINV-GFP and BQCV) and Apis mellifera rpl8 analyses consisted of a single pre-incubation 95ºC (3 min), 40 cycles of 95ºC (5 s), and 60ºC (20 s); primers listed in Supplementary Table S1. Positive and negative controls, including the use of RNA templates from no RT enzyme cDNA reactions, were included for all qPCR analyses and exhibited the expected results.
To quantify the viral RNA (i.e., genome and transcript) abundance in each sample target SINV-GFP qPCR amplicons were cloned into the pGEM-T (Promega) vector, as described in Flenniken and Andino et al. 2013 1 . Plasmid standards, containing from 10 9 to 10 3 copies per reaction, were used as qPCR templates to assess primer efficiency and generate the SINV-specific standard curve used to quantify the viral RNA copy numbers within this range of detection 1 . The qPCR primers for RNAseq validation were designed using Primer3Plus and typically designed to have 60ºC annealing temperatures 12 (Supplementary Table S1). Melt point analysis and 2% agarose gel electrophoresis ensured qPCR specificity 13 . Primer efficiencies were evaluated using qPCR assays of cDNA and plasmid dilution series, and calculated by plotting log 10 of the concentration versus the crossing point threshold (C(t)) values and using the primer efficiency equation, (10 (1/Slope)-1) x 100 (Supplementary Table S16).
The ΔΔC(t) method was used to calculate relative abundance of Sindbis-GFP in individual bees because it was most accurate since this method ensures that results are not skewed by inadvertent differences in RNA reverse-transcription efficiencies and starting cDNA template abundance 11,13,14 . Specifically, the ΔC(t) for each sample was calculated by subtracting the Am rpl8 C(t) from the SINV-GFP C(t). The honey bee gene encoding ribosomal protein 8, Am rpl8, was selected as an appropriate housekeeping gene for qPCR, since it has been utilized in several other studies [15][16][17][18] . Analysis of the RNASeq data presented herein confirmed that rpl8 expression was not significantly different in all sequenced libraries. The ΔΔC(t) was calculated by subtracting the average virus-infected ΔC(t) values from the ΔCt values for each treatment group. For host gene expression analyses and RNAseq validation, the percent gene expression for each gene of interest (GOI) was calculated using the following formula: 2^-ΔΔC(t) x100 = % gene expression, in which ΔC(t)=GOI C(t)-rpl8 C(t), and ΔΔC(t)= sample ΔC(t)mock-infected control ΔC(t).
To statistically examine the differences between the relative virus abundance of each treatment group for each time point, the SINV-GFP copy numbers and calculated relative abundance of each sample and was imported into R 19 . Based on previous work 1, 20 , we hypothesized that bees (n=10) co-injected with dsRNA or poly(I:C) would have decreased relative virus abundance as compared to the virus-only treated group. To examine relative virus abundance between treatment groups (e.g., virus-infected bees an dsRNA or poly(I:C) co-treated bees) that had equal variance and normal distribution we performed one-tailed students t-tests. Analysis of honey bee host gene expression revealed unequal variance between treatments groups and thus Welch's ttests were used to identify statistical differences in host gene expression.

RNAseq Library Preparation and Sequencing
Bees were obtained from honey bee colonies which are subject to naturally occurring infections that may confound transcriptional results 21 , so individual bee cDNA was screened for pre-existing infections via PCR for several honey bee pathogens (Supplementary Table S3) using the following PCR thermocycler protocol: 95 °C (5 min); 35 cycles of 95 °C (30 s), 57 °C (30 s), and 72 °C (30 s), followed by final elongation at 72 °C for 4 minutes. If the sample was positive for a pathogen, the quantity was then assessed using qPCR. The RNA isolated from the abdomens of at least three representative bees with low (< 2,000 DWV and/or BQCV virus genome copies versus 7 x 10 4 -7 x 10 6 SINV-GFP copies) to no pre-existing infections for each treatment group and time point were selected for transcriptome sequencing for a total of 47 individual bees (Supplementary Table S3). The abdomen was selected as the tissue of interest for RNASeq analysis of honey bee immune responses since virus infections would have naturally disseminated into this site post inoculation via intra-thoracic injection and it houses the fat body, which is an important tissue for immune function.
Prior to RNASeq library preparation, RNA from each sample was further purified using Qiagen RNeasy columns, including on-column DNase Treatment (Qiagen  Table S2), which is in the range of coverage reported in other honey bee transcriptome studies 3,4,24 . Previous work has shown that sequencing as few as 1 million reads can provide significant and valuable information of differential expression patterns in honey bees 24 . Sequence data was deposited into the NCBI Sequence Read Archive under accession number SRP101337 and is also linked with NCBI BioProject #PRJNA377749.
Specific code used for RNAseq analysis is listed below. Briefly, the programs FastQC and fastx-toolkit were used to remove low quality reads (<Q30 To further investigate the function of the DEGs, representative protein sequences (the longest sequence if there were splice variants) of every known honey bee gene were queried against the D. melanogaster protein database via reciprocal BLAST+ 31 to identify fruit fly orthologs and homologs of the honey bee genes there is a greater amount of gene ontology information for D. melanogaster genes compared to Apis mellifera genes. The honey bee genome encodes approximately 15,000 genes of which 13,592 genes are mapped and provided in the Amel4.5 genome annotation 23 . We annotated 8,944 genes (~66%) as homologs (of which 7,006 were reciprocal best hits or orthologs) to genes encoded by the fruit fly D. melanogaster genome, which encodes ~13,600 genes 23,32 . Biological processes (BP) functional enrichment analysis was performed with DAVID 33 .

Comparative analysis of DEGs in virus-infected bees
To identify the shared and unique DEGs in virus-infected bees, we compared our dataset to other studies that examined gene expression in virus-infected bees. We compared the genes that were differentially expressed in SINV-GFP infected bees 72 hpi to DEGs in symptomatic IAPV-fed bees 3 , SBV and DWV-infected bees 4 , adult honey bees naturally infected with IAPV 34 , and a common DEG list that was compiled from 19 gene expression data sets including Varroa destructor-parasitized and virus-infected bees 35 . We used NCBI gene ID as a common identifier because DEG lists were generated using different technologies (i.e., microarray and high throughput sequencing) and versions of the Apis mellifera genome and transcriptome. The DEGs in SBV and DWV-infected bees described by Ryabov et al., 4 were correlated to our reference gene list using the BeeBase Bee IDs, and resulted in the common identification of 1,315 of the 1,638 DEGs. Similarly, we assigned common gene identifiers for 662 of the 753 differentially regulated transcripts in the symptomatic IAPVfed bees 3 . The DEGs identified in naturally IAPV infected adult honey bees using microarray analysis 34 was kindly provided by Dr. Judy Chen. This list consisted of microarray probe IDs, therefore, in order to compare the microarray data to our data set, the microarray oligonucleotide probe sequences were aligned to Apis mellifera refseq RNAs using BLAST (megablast, plus / plus strands) and the corresponding NCBI Entrez Gene IDs were identified. We identified 4,402 NCBI Entrez Gene IDs corresponding to the 5,637 probes for specific transcripts and transcript variants. Of the 4,402 probes for which we identified 4,307 Gene IDs, 3,109 of them matched distinct genes that were also in our gene list. The NCBI Entrez Gene IDs of each DEG list (Supplementary Table  S14) were then compared via Venn diagram analysis 36 (Supplementary Table S15). Pairwise comparisons between studies identified shared DEGs and the statistical significance of gene overlap was assessed using hypergeometric tests 37 (Supplementary Table S15).

Identification of previously unrecognized transcript
RNASeq analysis determined that reads aligning to LOC725387 were more abundant in virus-infected bees. To identify the gene or genes encoded by these differentially expressed reads, the consensus nucleotide sequence was used to query the NCBI Nucleotide collection (nr/nt) and A. mellifera databases using blastn 38 , Sanger sequencing was performed to verify transcript sequence and length, and the results were evaluated using Geneious 39 . Together, these analyses revealed that we identified a previously unrecognized transcript, A. mellifera probable cyclin-dependent serine/threonine-protein kinase (MF116383, 5,158 nt), which is longer the originally annotated A. mellifera probable serine/threonine-protein kinase clkA (LOC725387, XM_001121241.4, 1,403 nt).
In brief, we utilized LOC725387 RNASeq consensus sequence to query the NCBI Nucleotide nr/nt data base and identified an A. cerana transcript annotated as a probable cyclin-dependent serine/threonine-protein kinase DDB_G0292550 (LOC107994302, XM_017051141.1) as the top blastn result, which contained 95% of the submitted sequence and shared 91% identity (E-value = 0, 95% query coverage, 91% identity, 1-6% gaps); additional top blastn hits included A. dorsata GATA zinc finger domain-containing protein 14-like. When the LOC725387 RNASeq consensus sequence was used to query the A. mellifera database, the top blastn result only covered 24% of the query sequence (i.e., A. mellifera probable serine/threonine-protein kinase clkA, XM_001121241.4; E-value = 0, 24% query coverage, 99% identity, 0% gaps). To further characterize the LOC725387 transcript, we Sanger sequenced 5,027 nts (~2-3x coverage) and obtained the most 5'end of this transcript from RNASeq data (131 bp, >2,000x coverage) (Supplementary Table S1 and Fig. S9). Together nucleotide and amino acid alignments indicate the RNAseq reads aligning to LOC725387 are most similar to a computationally predicted A. cerana cyclin-dependent serine/threonineprotein kinase DDB_G0292550 (Supplementary Figure S9), therefore, we refer to the gene identified herein as A. mellifera probable cyclin-dependent serine/threonineprotein kinase (Supplementary Fig. S9) and submitted the sequence of this transcript to NCBI (MF116383).

Vennt
a. Download python script "venn.py" from http://drpowell.github.io/vennt/ b. Works with Cuffdiff output and generates list of DEGs based on manually set q values, for which it was set to <0.05.