Cerana

Raw sequencing reads with phred scores ≤ 20 were filtered out using the CLC_quality_trim (CLC 3.22.55705). Duplicate sequences were removed with the remove_duplicate program (CLC-bio) using the default options. After filtration, genome libraries with inserts of 500 bp, 3 kb, and 10 kb were assembled using the AllPaths-LG (version 42411, [31 (link)]) algorithm with default parameters. The A. cerana genome sequence is available from the NCBI with project accession PRJNA235974. Repeat elements in the A. cerana genome were identified using RepeatModeler (version 1.0.7, [98 (link)]) with default options. Subsequently, RepeatMasker (version 4.03, [99 (link)]) was used to screen DNA sequences against RepBase (update 20130422, [100 (link)]), the repeat database, and mask all regions that matched known repetitive elements. Comparison of experimental mitochondrial DNA to published mitochondrial DNA (NCBI accession GQ162109) was performed using the CGView Server with the default options [101 (link)]. The percent identity shared between the A. cerana mitochondrial genome assembly and NCBI GQ162109 was determined by BLAST2 [102 (link)]. To examine the distribution of observed to expected (o/e) CpG ratios in protein coding sequences of A. cerana, we used in-house perl scripts to calculate normalized CpG o/e values [43 (link)]. Normalized CpG was calculated using the formula:

where freq(CpG) is the frequency of CpG, freq(C) is the frequency of C and freq(G) is the frequency of G observed in a CDS sequence.

To compare tissue transcriptomes of A. cerana and A. mellifera, RNAseq reads from antenna tissue of each species were mapped against annotated Gr, Or, and Ir gene sequences using RSEM [113 (link)]. TMM (Trimmed Mean of M-values)-normalized FPKM (fragments per kilobase per million) expression values were calculated by running perl scripts included in Trinity assembler (version r2013-02-25) and were used to draw a heat map [98 (link)].

Microsynteny analyses of orthologous Or genes between A. cerana and A. mellifera were carried out based on a phylogenetic tree. In addition, we performed reciprocal BLASTZ searches [114 (link)] between two species in syntenic regions and identified clear orthologs based on E-values (cut-off: 1e⁻⁴⁰) and phylogenetic analyses.

The nucleotide and amino acid sequences for tobacco (Nicotiana tabacum) KED [6 (link)] were used as a starting point for the initial search using BLASTN and BLASTP tools against the GenBank and OneKP database including non-redundant nucleotide and protein sequences, whole-genome shot gun, expressed sequence tags, high throughput genomic sequences, UniProtKB, transcriptome shotgun assembly proteins and protein data bank. Initial search using both the coding nucleotide sequences and the amino acid sequences identified 32 eudicots and at least one monocot (Elaeis guineensis). Subsequent searches were performed against the E. guineensis amino acid sequence through Liliopsida (monocotyledons) database to identify matching sequences of monocots. Likewise, using retrieved sequences to systematically search databases of the same orders and families of eudicotyledons yielded more species of possible matching sequences. Similar strategy was used to identified KED sequences from gymnosperms. Further searches were done sequentially by narrowing organism groups to find matches from more closely related species. However, it must be pointed out that the database search was aimed at surveying broadly the possible taxonomic presence of the KED gene and the retrieved sequences are not by no means an exhaustive outcome due to genomic sequence availability and the annotation quality of the public databases.
To search for possible KED-rich sequences in charophytes, bacteria and animals, KED protein and nucleotide sequences from plants were first repeatedly blasted through each of the intended organism groups in the databases. Then each match was further examined by retrieving its sequence from the database. Translation tool was used to generate open reading frames, followed by amino acid composition analysis, specifically for K+E+D content, to score the putative KED candidates. Once a KED sequence was identified from one taxon group (for example, charophyte), this sequence was used to search the entire available entries from this group. This way, sequences predicting KED-rich open reading frames in genomes of several charophyte, bacterial and animal species were identified.
During the course of searching animal KED candidates, a 6,229-amino acid microtubule-associated protein futsch from honeybee (Apis cerana) was found to contain an internal KED-rich region, whereas its N- and C-terminus portions have normal K, E and D contents. To illustrate examples of the presence of KED sequences in animal species, this 750-amino acid internal KED-rich region was arbitrarily taken out for demonstration in this study.
All retrieved sequences of possible matches were manually reviewed and verified for proper open reading frames and translated sequences. Wherever applicable, both genomic sequences and mRNA sequences were matched to verify the correct coding sequences. The full-length, translated sequences with considerable sequence identity and a high percentage of KED (K+E+D% greater than 30%) were designated as a candidate match.
Only partial KED sequences were available for two plants: cedar (Cryptomeria japonica, a gymnosperm; without C-terminus) and barley (Hordeum vulgare, a monocot, angiosperm; without N-terminus). However, they both still possessed the conserved domain (see “Results” below), therefore were included in sequence comparison analysis. But because their KED protein lengths were unknown and would distort the analysis parameters, they were excluded from the dataset for phylogenetic analysis described below.

SMC +  + 1.13.1 [61 (link)] was used to estimate the changes in Ne over the past one million years (Ma). The mutation rate was set to 5.27 × 10⁻⁹ per base pair per generation, following a divergence estimate of 7 Ma between A. mellifera and A. cerana [62 (link)]. The generation time was assumed to be one year. The polarization error was set to 0.5.

Ten adult A. cerana workers from each colony were pooled and homogenized as described for the initial RNA extraction step above. Total genomic DNA was extracted using a DNA purification kit (PureLink Genomic DNA Mini Kit, Invitrogen, Carlsband, CA, USA) according to the manufacturer’s instructions. DNA samples were stored at −20 °C prior to molecular screening for microsporidia, fungi, and bacteria.