Differentially and Co-expressed Genes in Embryo, Germ-Line and Somatic Tissues of Tribolium castaneum

Transcriptomic studies of Tribolium castaneum have led to significant advances in our understanding of co-regulation and differential expression of genes in development. However, previously used microarray approaches have covered only a subset of known genes. The aim of this study was to investigate gene expression patterns of beetle embryo, germ-line and somatic tissues. We identified 12,302 expressed genes and determined differentially expressed up and down-regulated genes among all samples. For example, 1624 and 3639 genes were differentially increased in expression greater than or equal to twofold change (FDR < 0.01) in testis vs. ovary (virgin female) and ovary vs. embryo (0-5 hr), respectively. Of these, many developmental, somatic and germ-line differentially expressed genes were identified. Furthermore, many maternally deposited transcripts were identified, whose expression either decreased rapidly or persisted during embryogenesis. Genes with the largest change in expression were predominantly decreased during early embryogenesis as compared to ovary or were increased in testis compared to embryo. We also identify zygotic genes induced after fertilization. The genome wide variation in transcript regulation in maternal and zygotic genes could provide additional information on how the anterior posterior axis formation is established in Tribolium embryos as compared to Drosophila. Together, our data will facilitate studies of comparative developmental biology as well as help identify candidate genes for identifying cis-elements to drive transgenic constructs.

differences between males and females in insects has been performed but is variable (Papathanos et al. 2009). Some sexually dimorphic characteristics, such as body size or development rate are also influenced by natural variation, thus regular adjustment and recalibration are required for such systems to be used. For the early sterile insect technique (SIT) programs for insects, especially Aedes aegypti, sexes were separated using differences in pupal body size. However, many of these methods are not directly transferable to other insect species.
As an alternative, transgenic sexing systems (transgenic animals show male or female specific lethality during embryonic and early larval stages, leading to either a male or female specific population) can be developed using numerous technological approaches for the purpose of efficient sex separation in insects. Unlike biological methods for sexing, genetic methods may be more broadly applicable across species; with some minor modifications, the constructs developed for one species can be adopted to another unrelated species. For example, modifications to specific individual components of the construct such as endogenous species-specific cis elements or sex-specific splicing variants and lethal effectors may increase the efficiency of the sexing construct (to develop a construct in which a transgene leads to a high level of expression in either male or female), although simultaneously avoiding the potential negative effects of accidental species transfer. For insect transformation three components are important, (1) a gene vector (transformation method), (2) an effector gene and (3) cis elements, particularly the promoter/enhancer sequences which control transcription of the gene in a specific tissue at specific a developmental time point such as sex-specific expression of a lethal gene leading to elimination of one sex and population suppression. However, genetic sexing in a targeted species requires known differentially expressed genes in the different sexes.
Coleoptera (Beetles) are the most diverse animal group on earth and contain one fourth of all species described and includes many major pests of crop plants (Hunt et al. 2007). T. castaneum is a genetically tractable model beetle species with a published whole genome sequence (Tribolium Genome Sequencing Consortium et al. 2008). Tribolium follows the XX/XY sex determination system (Male XY and Female XX), and male and female Tribolium have some sexually dimorphic characters such as black spots on the first pair of legs of male adults (male beetle sexual dimorphism) which are absent in females, as well as differences in the appearance of male and female pupae (beetle pupal sexual dimorphism). To screen and separate beetles on a large scale based on these naturally occurring biological phenotypes is likely not feasible and certainly very difficult. Therefore, a study of gonadal differentiation and embryo development gene expression in Tribolium would be useful for the development of sexing strategies for this and other coleopterans insects important in agriculture and medicine.
Deep sequencing of mRNA (RNAseq) has been successfully used for differential gene expression analysis in a wide variety of species and conditions (Akbari et al. 2013;Graveley et al. 2011;Guo et al. 2018;Ruan et al. 2018). We performed RNAseq analysis of transcripts isolated from testis, ovary, carcasses and early embryos in Tribolium castaneum. We detected expression from a total of 12,302 genes in one or more of these tissues, corresponding to 75% of the total annotated genes (16,593) in the Tribolium genome. The identification of differentially or unique expressed genes reported here will facilitate future work to identify sex and early embryo specific cis elements for heterologous gene expression and to improve genetic sexing in Tribolium allowing for embryogenesis studies and improved pest control in the field.

MATERIALS AND METHODS
Red beetle strain and culture The T. castaneum strain was obtained from Jeff Demuth at the University of Texas at Arlington. Beetles were reared on flour medium (95% flour, 5% yeast by weight), and caged in glass jars with tight-fitting fine mesh closures. Beetles were housed in a growth chamber at 25°w ith 60-80% relative humidity and 12/12 hr light cycling. Populations of beetles were moved to fresh flour medium once per month with initial population densities of approximately 122 beetles/1 g flour medium.
RNA isolation and RNAseq library preparation Beetle pupae were separated based on sex prior to eclosion. Adult beetles (1-2 weeks post-eclosion) were used for embryo collection. Adult beetles were placed on flour pre-sifted with no. 50 sieve to lay eggs at 25°. The adults were removed by sifting with no. 50 sieve and transferred to fresh flour to continue embryo collection. The embryos were collected from the flour by sifting with a no. 30 sieve. Embryos were collected at 0-5 hr and 6-11 hr, respectively, in three independent biological replicates (each replicate had 80-100 embryos and were processed in a single day per sample). Similarly, ovary, testis and carcasses from male and female adult beetles at four days post eclosion (30 beetles in each replicate and processed in a single day per sample) were collected to evaluate transcriptional activity. Samples were snap frozen in liquid nitrogen and then transferred to -80°prior to RNA extraction. RNAextraction and library preparation were done simultaneously with all samples. Total RNA from these tissues was extracted with Trizol (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. After total RNA extraction, RNA was treated with Turbo DNase (Ambion/Applied Biosystems, Austin, TX). RNAseq libraries were prepared using NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB #E7490) for mRNA isolation and NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #E7530S) according to the manufacturing instructions. Each library had a peak size of approximately 300 base pairs. cDNA library quality was assessed on an Agilent 2200 TapeStation using a High Sensitivity D1000 tape. Libraries were quantified using the Invitrogen Qubit 2.0 Fluorometer High Sensitivity dsDNA assay. All libraries were normalized to a 4 nM concentration and then pooled in equal volumes for denaturing and diluting for sequencing. The library pool was denatured and diluted to a final concentration of 1.8 pM following the Illumina NextSeq Denature and Dilution protocol including a 1% PhiX DNA spike in. Sequencing was performed on an Illumina NextSeq 500, 75 cycles, high-output sequencing run generating approximately 550 Million sequencing reads. We used Illumina genome sequencing platform (NextSeq500 illumina R ) at the Texas A&M Institute for Genome Sciences and Society Shared Molecular Genomics Core to sequence the aforementioned samples, producing unpaired 75-nt raw reads.

Data Analysis
Sequencing bcl files were automatically uploaded to Illumina's cloud service, BaseSpace. BaseSpace automatically generates FASTQ files, demultiplexes and trims adapter sequences. FastQC (version 0.11.5) tool was used to check the quality of the reads in each sample. HiSat2 software version 2.0.5 (Kim et al. 2015) was used to align reads to the Tribolium genome (GCA_000002335.3_Tcas5.2_genomic_RefSe-qIDs.gff, obtained from www.ncbi.nlm.nih.gov) using default parameters. BEDTools (version 2.27.0) was used to produce raw counts. We used two approaches to determine which genes encode maternal, zygotic, soma or germ-line unique or increased gene expression. First, we identified genes whose expression is co-regulated across samples by using soft clustering (50% mem Soft Cluster graph) [Mfuzz package (Kumar and Futschik 2007) in 'R' program]. The "cselection" function was used in Mfuzz to determine cluster number. We used the method proposed by Schwaemmle and Jensen (Schwämmle and Jensen 2010) to estimate the value of fuzziness parameter "m", which is based on the distribution of distances between genes in a given data set, to set a suitable value for "m". Maternal genes were classified based on gene expression in both the Tribolium ovary and early embryos. Zygotic transcripts and induced zygotic transcripts were identified as mRNAs that were co-regulated in embryo and did not appear in other tested samples. Second, we selected the set of genes whose expression changes greater than or equal to twofold in tissue or embryo in pairwise comparison to other tested samples (edgeR in R program, FDR , 0.01). Edge R software (McCarthy et al. 2012) in "R" package was used to normalized the libraries using the trimmed mean of M-values method and for differential gene expression analysis. Genes were kept in the analysis only if the counts per million (CPM) was equal or greater than five in three or more samples.

Quantitative real-time PCR
The software Primer-3 (http://frodo.wi.mit.edu/) was used to design primers for qPCR analysis. Primers were designed by the rules of highest maximum efficiency, and sensitivity rules were followed to avoid formation of self and hetero-dimers, hairpins and self-complementarity. The primer sequences used in this study are shown in Table S1. In brief, single-strand cDNA was synthesized as follows: 500 ng of total RNA in 11 ml of sterile water previously treated with DNase using the DNA-free kit (Ambion www.ambion.com) following the manufacturer's instructions. The reverse-transcriptase reaction to generate the cDNA for use in quantitative real-time PCR was carried out using the First Strand cDNA Synthesis kit (Fermentas) as follows: 1 ml of oligo d (T) primer was added to the 11 ml of total RNA. The mixture was heated at 65°for 5 min, and then placed on ice, and the following were added: 4 ml of 5X first-strand buffer, 2 ml of dNTPs, 1 ml of RNase inhibitor, and 1 ml of reverse transcriptase. cDNA synthesis was performed at 42°for 30 min and 50°for 60 min. Reactions were stopped by heating samples at 95°for 2 min.
qRT-PCR was performed in optical 96-well plates on a BioRad Real-Time PCR Detection System using the Absolute QPCR SYBR green Mix (ABgene) to monitor double-stranded DNA synthesis in combination with ROX as a passive reference dye included in the PCR master mix. Amplification conditions were 10 min at 95°to activate the polymerase, followed by 40 cycles at 95°for 30s, 60°for 30s and 72°for 30s. A melting curve analysis was performed in order to verify that the amplicon was at a correct T m and therefore of the correct length of the predicted transcript. Results were normalized to the mRNA level of T. castaneum actin (ACT) and ribosomal protein RpSL32 as housekeeping genes (Table S1), and data calculated according to the deltadelta Ct method. Each experiment was repeated with three independent mRNA samples (biological replicates), and each reaction was repeated three times to minimize intra-experiment variation (technical replicates). All of the results were analyzed with CFX manager software.

Drosophila melanogaster homologs in T. castaneum
The iBeetle web server (http://ibeetle-base.uni-goettingen.de/) (Dönitz et al. 2018) was used to identify T. castaneum homologs to the D. melanogaster genes, which calculated homologs based on OrthoDB. Lists of D. melanogaster known maternal and zygotic genes, also genes highly expressed in the D. melanogaster germ cells and somatic tissue were retrieved from Table S12-30 of (Arbeitman et al. 2002), Table S4, S6, S7 and S8 of (De Renzis et al. 2007), additional files 3, 8 and 9 of (Lebo et al. 2009) and additional data 21-23 of (Parisi et al. 2004). Additional lists of T. castaneum genes expressed greater than or equal to twofold change in fertilized and unfertilized eggs were obtained from additional file 1 of (Preuss et al. 2012). In this study, we selected the D. melanogaster homologs of genes expressed in the Tribolium embryo and different tissues, which either showed unique gene expression pattern in clusters by using soft clustering analysis or in pairwise comparisons showed greater than or equal to twofold change in gene expression (FDR ,0.01) in edgeR analysis in respective tissues.
Gene Ontology (GO) and functional annotation Functional activity annotation was performed on genes in soft clusters and differentially expressed genes greater than or equal to twofold change by using the g:Profiler (g:GOSt Gene Group Functional Profiling) Gene Ontology web server (Reimand et al. 2016). To extract GO terms information and enrichment analysis for the selected genes, we used the following parameters, user threshold (p-value) = 0.001. Significance threshold: g: SCS threshold and organism panel was selected as T. castaneum.

Data availability
All raw and processed RNA seq data from the study are available at the gene expression omnibus (GEO) (Edgar et al. 2002)

RESULTS AND DISCUSSION
Here, we report transcriptional profiles during embryo development, male carcass (lacking testis), female carcass (lacking ovary), testis and ovary by deep RNA sequence analysis (Wang et al. 2009). We obtained transcriptome reads for 0-5 hr embryos (average from three replicates 33,458,920 reads); 6-11 hr embryos (average from three replicates 29,034,300 reads); male testis (average from two replicates 29,359,190 reads); male carcass (average from three replicates 32,279,190 reads); female ovary (average from three replicates 31,321,690 reads) and female carcass (average from three replicates 32,059,370 reads) as shown in Figure S1. Out of 16,363 total annotated genes, 12,302 were found to have a normalized expression (CPM) greater than or equal to 5 in three or more samples. A Multidimensional Scaling Plot (MDS) of these samples using edgeR (Robinson and Smyth 2008) on the normalized counts showed similarity within sample sets and dissimilarity of gene expression among sample sets, indicating consistency among replicate samples ( Figure S2).
To examine the global dynamics of differential gene expression in the Tribolium embryo, germ-line and somatic tissues, we compared the normalized expression levels for every sample with every other sample resulting in a global view of gene expression as shown via clustering map (Figure 1). The clustering pattern observed showed that gene expression in male and female germ tissues are substantially different from each other and to both somatic tissues (carcass lacking ovary or testis) and the early embryo. In contrast, the transcriptional profiles of carcass tissues isolated from males and females were highly similar to each other, suggesting minimal sexually dimorphic expression in the adult soma. Similarly, the expression profiles of two embryo stages were similar to each other.
To better understand the major patterns of gene expression, we performed a differential expression analysis between all pairs of samples in the embryo and among tested tissues, classifying individual genes as differentially expressed if they exhibited a twofold change at a false discovery rate of ,0.01 (FDR ,0.01) as given in Table S2. Additionally, we used a soft clustering algorithm [Mfuzz package (Kumar and Futschik 2007) in 'R' program] and established 15 groups of co-regulated genes, with each cluster including between 312 and 1744 expressed genes ( Figure 2, Table S3). In the following sections, we present the results from these analyses focusing on each developmental stage/tissue analyzed.

Identification of maternal gene expression in the Tribolium early embryo
We identified genes differentially expressed in the embryo (0-5 hr) as compared to the ovary, testis, embryo (6-11 hr), and male and female carcasses (Table S2). In pairwise comparisons of embryo (0-5 hr) to other tested samples as given in Table S4, 35 genes were found to be preferentially expressed in the early embryo (0-5 hr). Out of these 35 genes, twelve had homologs of maternal transcripts in the fly embryo (Table S4), including the TC032769 (caudal) a homolog of the fly gene (FBgn0000251), which is highly expressed in the fly embryo and involved in anterior/posterior patterning formation in the fly embryo. When compared only to the ovary, 62 genes were expressed more than twofold higher in embryo (0-5 hr). Out of these, 20 genes had fly homologs that were also expressed in the fly embryo (Table S4). In strong contrast, comparing gene expression between the early embryo (0-5 hr) to ovary revealed 3639 genes with reduced expression (Table  S5). Of these, 1833 genes had a more than fourfold decrease expression in the early embryo (0-5 hr) as compared to ovary, with 27 genes decreased by more than tenfold (Table S5). These genes largely correspond to Cluster-2, which included 1331 genes highly expressed in the ovary (Figure 2), and include those important for the formation and maintenance of this tissue as well as maternally synthesized products required for oocyte formation and early embryogenesis.
As the best studied model system for insect development outside of Drosophila, we sought to better identify homologs of maternally deposited Tribolium genes that were also maternal in Drosophila (Preuss et al. 2012;Arbeitman et al. 2002) (Figure 3). Among Clusters-2, -7, and -10 we identified 3,817 Tribolium genes with expression patterns consistent with those of maternal genes, or 31.02% (maternally deposited) of total transcripts sequenced in this study. In particular, Cluster-7 was assigned genes with higher expression in the Tribolium early embryo (0-5 hr) and ovary, but decreased expression after 5 hr (Figure 2). For comparison in Drosophila, most maternal genes decrease expression by more than threefold within 6.5 hr of embryonic development or zygotic expression initiation (Arbeitman et al. 2002). From prior published studies, we identified 1,439 known maternal genes from Drosophila and used these to identify homologs in our dataset, specifically in Cluster-2, 7 and Cluster-10 ( Figure 3). Out of 1,439 known Drosophila maternal genes, 71 and 123 homologs of known Drosophila maternal genes were identified in Cluster-2 and -7 respectively (Figure 3). For example, TC002055 present in Cluster-7 is a homolog of D. melanogaster FBgn0026181, which is expressed during oogenesis in ovary and deposited to the egg (Verdier et al. 2006). TC015003, is a homolog of the Drosophila FBgn0026206 (meiotic P26) (GO:0051321, meiotic cell cycle), which is expressed during oogenesis and was classified as a maternally loaded expressed gene in the fly embryo (Arbeitman et al. 2002). Similarly, the Tribolium intersex gene (sex determination) had high expression in early embryo (Figure 2, Cluster-7), and the fly homolog of this gene, FBgn0001276, is also highly expressed in the fly embryo (Arbeitman et al. 2002). Cluster-10 was assigned 1006 genes expressed in both embryo stages and also in the ovary. Of these, 89 were found to have homologs to known maternal genes that persisted during embryogenesis (Figure 3, Cluster-10). For example, myd88 (FBgn0033402) is a known fly maternal gene present in Cluster-10, myd88 (TC003185) is co-expressed in the Tribolium early embryo, zygotic stage and ovary.
Differential expression of genes in the Tribolium zygote While comparisons between transcript levels in ovary and early embryos can reveal maternally-derived gene products, comparisons with later stage embryos were designed to reveal early zygotic transcripts. At 25°, the embryonic development time of T. castaneum is significantly longer than that of Drosophila, and the Tribolium embryo at less than 5 hr is referred as pre-zygotic (Ninova et al. 2016). When we compared postzygotic embryo (6-11 hr) and pre-zygotic embryo (0-5 hr) samples, 719 genes were identified as early zygotic genes in the Tribolium embryos (6-11 hr), with increased expression as compared to early embryo (0-5 hr) (Table S6; Figure S3). This includes gene TC005375, whose transcripts were substantially up-regulated in embryo (6-11 hr) compared to all tested samples (Table S6). TC005375 is a homolog of D. melanogaster FBgn0002565 (Larval serum protein2) in Tribolium. Additionally, expression of TC015380 (homolog of D. melanogaster FBgn0033464) increased more than fourfold in 6-11 hr embryos in Figure 1 Gene expression pattern revealed by clustering analysis in early embryo and different tissues. Clustering of normalized differentially expressed genes in Tribolium based on heat map analysis made using pheatmap package in 'R'. Euclidean distance was used to estimate pairwise sample distances and sample clustering was done using the k-means algorithm. The columns of the heat map figure indicate samples and the rows represent different genes. The bar color reflects the gene expression levels as TMM normalized log-CPM. Color key indicates the intensity associated with normalized expression values. Red and Blue colors indicate higher expression and lower expression respectively. comparison 0-5 hr embryos. In the Tribolium embryo (6-11 hr), fifteen genes were identified with significantly increased expression in pairwise comparisons with each of the other tested samples (Table S6); out of these, nine had homologs that are early zygotic genes in the fly (Table S6). Genes that are expressed constitutively, particularly in the early embryo, could be potential candidates to use for gene drive systems in insects (Buchman et al. 2018).
Cluster-10 (n = 1006 genes) included many genes with high expression in early embryo (0-5 hr), embryo (6-11 hr) and ovary samples. This may represent the persistence of maternal transcripts or, more likely, continued expression during the early zygotic phase.
We identified many genes with expression in both early embryos (0-5 hr) and embryos (6-11 hr) as shown in Figure 2, Cluster-8, which are not expressed in other tested samples. Known zygotic genes were used as queries to identify the any homologs of zygotic expressed genes in Tribolium. In cluster analysis, the expression profile pattern of 456 Tribolium genes in the embryo samples was similar to that of 137 known zygotic genes in D. melanogaster (Figure 2; Cluster-8; Table S6). Five previously characterized early zygotic genes (TC030997, TC008514, TC009052, TC015380 and TC033031) present in Cluster-8 were significantly up-regulated in the Tribolium embryos (6-11 hr).

Identification of genes highly expressed in the Tribolium male germ-line
Cluster-14 was assigned 531 genes, with expression observed primarily in the testis (Figure 2, Table S3). Pairwise comparisons with other samples ( Figure 4A-J) revealed 1,014 genes with increased expression in the testis when compared to the ovary, female and male carcasses ( Figure 5, Table S7). Genes highly expressed in the Tribolium testis were examined for homologs from D. melanogaster, leading to the identification of 272 genes with orthologs known to be expressed preferentially in the testis-of the fly (Table S7). Differential high expression of genes in the testis may play specific roles in spermatogenesis and development. For example, the gene TC032055 had the highest expression (greater than tenfold change) in testis as compared to the other samples as shown in Table S7. This gene is a homolog of D. melanogaster FBgn0036093 that has a domain of unknown function. FBgn0036093 has two transcripts and two different polypeptides and has expression observed during the pupal period and in adult males. The Tribolium homolog has two transcripts (XM_008203121.2 and XR_001575917.1).
Similarly, another gene highly expressed in the testis was TC006703, which codes for the DNA double-strand break repair Rad50 ATPase. Transcript abundance for this gene ranges from a peak of moderately high expression in testis to very low expression in female carcass (Table S7). Rad50 ATPase is an essential component of MRN Figure 2 Stage specific changes in gene expression. Genes co-expressed in the Tribolium different embryo stages and tested tissues. The X-axis represents samples, samples 1-3 (embryo 0-5 hr), 4-6 (embryo 6-11 hr), 7-9 (female carcass), 10-12 (male carcass), 13-15 (ovary), and 16-17 (testis), with the Y-axis representing change in expression. Line color was assigned to each gene, with red (1) indicating high association in membership to the cluster. Purple, green or blue lines correspond to genes with lower membership values; "n" indicates number of genes assigned to each cluster.

Figure 3
Cluster of highly co-expressed genes in the Tribolium ovary and embryo. Genes expressed and grouped together in the Tribolium ovary (maternal genes) is shown in Cluster-2 (yellow filled color oval). Genes co-expressed in the Tribolium embryo (0-5 hr) and ovary (maternal genes and reduced in expression rapidly) is shown in the Clusters-7 (green filled color oval). Genes co-expressed in both embryo stages and ovary is grouped together in Cluster-10, coral filled color oval (maternal genes). The purple filled color oval represents the known maternal genes. The overlap represents the known maternal genes in Tribolium.

Figure 4
Gene expression patterns in the Tribolium testis. Differential genes expression in the Tribolium testis compared to male carcass (A -B), female carcass (C-D) ovary (E-F), embryo (0-5 hr) (G-H) and embryo (6-11 hr) (I-J). Colored dots represent individual differentially expressed genes. Red and blue color dots indicate significantly differentially expressed genes; red colored dots on Volcano plots are significant and above twofold or below twofold. N= number of genes DE in greater than or equal to twofold change or lesser than or equal to twofold change.
(Mre11-Rad50-Nbs1) complex; its molecular function is involved in 39-59 exonuclease activity and in biological process double-strand break repair via homologous recombination. The MRN complex contributes to processes such as recombination, DNA damage repair, genome stability and meiosis in germ cells. In Drosophila, Rad50 also forms a complex with Mre11 to cap telomeres during embryogenesis. Its abundance in the Tribolium testis may signify a special need for homology-based DNA repair during spermatogenesis.
Eleven Tubulin family members from the genome of T. castaneum were previously described; these were characterized as alpha tubulins, beta tubulins, gamma tubulins, delta tubulins and epsilon tubulins (Siebert et al. 2008). b 2 -tubulin specific cis-elements from different insects such as D. melanogaster, B. mori, A. aegypti and A. gambiae have been used to drive transgene expression in the testis (Michiels et al. 1989;Xu et al. 2015;Smith et al. 2007;Catteruccia et al. 2005). We identified four different b-tubulin homologs TC009589, TC034766, TC010829, and TC009035 from T. castaneum reference genome ( Figure  S4A). The genetic structure of these four b-tubulins homologs from T. castaneum is described in the Figure S4B. We found that only TC009035 was highly expressed (ninefold over other tissues) in testis (Table S7). TC009035 is a homolog of D. melanogaster b-Tub85D (CG9359) and CG9222 (FBgn0031784), with the b-Tub85D gene highly expressed in the D. melanogaster testis (Graveley et al. 2011). In Drosophila, this gene is transcribed in late third larval instar before the onset of meiosis in the developing testis and remains active throughout adulthood (Kemphues et al. 1982). A lower level of expression for b-Tub85D was also reported in other tissues such as in adult carcass and larval fat body (Gelbart and Emmert 2013).
To validate RNAseq data for genes expressed predominantly in the testis, we performed qPCR analysis as an independent measure of gene expression levels for TcRad50, Tc-b 2 -tubulin and an enolase gene.
In each case, the qPCR data analysis demonstrated robust expression in the Tribolium testis (Figure 6 A-D), confirming our RNAseq data and analysis.
Differential gene expression in the Tribolium female germ-line As stated above, the analysis of genes assigned to Cluster-2 revealed a group of genes that are transcribed in ovary and deposited into embryo but rapidly reduced in expression in embryo. Cluster-2 was assigned 1331 genes, out of which 71 were Tribolium homologs of fly femalebiased genes, which were transcribed only in the ovary and may have a role in ovary tissue development and oocyte synthesis. Two genes out of these 71 genes were TC009818 (homolog of D. melanogaster FBgn0020278; loco) and TC013108 (homolog of D. melanogaster FBgn0003187; quail) which have a role in oocyte dorsal/ventral axis specification and ovarian nurse cell to oocyte transport respectively. D. melanogaster loco may be required for nurse cell dumping during oogenesis (Pathirana et al. 2001). Mutations in D. melanogaster quail result in female sterility (Mahajan-Miklos and Cooley 1994). Our data suggests that the homologs of these two genes may be needed for egg production during oogenesis (Figure 2, Cluster-2). As an alternative method to identify differentially expressed genes in ovary, we compared gene expression profiles in ovary to all tested samples ( Figure  S5A-H), considering a threshold of greater than or equal to twofold change in gene expression in ovary as compared to testis, female and male carcasses at an FDR value , 0.01 (Table S8). By this method, 281 genes were identified as differentially expressed in the ovary tissue as compared to all other samples (Figure 7; Table S8). Such genes are good candidates for donating cis-acting control elements. For example, TC034053, which has greater than threefold (q-value 8.99x10 218 ) increased expression in the Tribolium ovary compared to testis, is a homolog of the D. melanogaster gene FBgn0039972 (meiosis arrest female protein 1), which has known expression in adult female and early embryo (Table S8).
Identification of genes differentially expressed in the Tribolium germ-line of either sex Sexual reproduction involves specialized variation of certain cellular processes in germ-line tissues such as sperm and egg production. To identify genes differentially regulated in the Tribolium germ-line, we compared expression levels of each expressed gene in the male germline tissues (testis) to expression in both male carcass and female carcass; simultaneously we compared gene expression levels in the ovary to expression levels in both male and females carcasses ( Figure 8A). From these comparisons, 374 genes had higher expression in the germ-line greater than or equal to twofold in both sexes relative to somatic tissues ( Figure 8A; Table S9). Of these germ-line expressed genes, 81 had homologs in Drosophila that are also preferentially expressed in the germ-line (Table S9). For example, TC001383 (fourfold change; q-value 2.91x 10 223 ) is a homolog of a D. melanogaster gene (FBgn0001086), which is involved in the anaphase promoting complex in mitosis that functions to regulate the three mitotic cyclins A, B and B3 in the egg and drive anaphase progression during meiosis (Frise et al. 2010).
From the same comparisons, we also sought to identify genes with significantly decreased expression in germ-line tissues as compared to the rest of the body ( Figure 8B, Table S10). We identified 939 genes, with significantly less expression in the male and female germ-line as compared to the dissected carcasses without these organs; with 27 genes more highly expressed only in the male carcass and 29 genes expressed preferentially in the female carcass ( Figure 8B and Table S10). Genes expressed preferentially in Tribolium somatic tissues were compared to their D. melanogaster homologs, leading to the identification of 264 homologs from D. melanogaster with corresponding soma-biased genes in Tribolium (Table S10). Among these, the expressions of 10 D. melanogaster homologs were also male-biased and four D. melanogaster homologs were female biased (Table S10). For example, TC000087 is a homolog of the FBgn0001208 (henna), which has a role in the biosynthesis of pteridin eye pigment. Similarly, TC010476 (greater than sixfold change; q-value 4.59X10 214 ) is a homolog of the D. melanogaster gene FBgn0005633 (flightin). Flightin (Fln) a flight muscle specific protein in D. melanogaster, and is not expressed in other muscle types and is involved for the correct assembly of the myosin thick filament. Clusters-1 and -15 include genes with expression in adult somatic tissues, while Cluster-3 (313 genes) had the smallest number of genes of all the clusters (Figure 2, Table S3) and was assigned genes with high expression in male carcass (lacking testis). Finally, Cluster-12 included 1057 genes ( Figure 2, Table S3), which were highly expressed in tissues from the adult stages analyzed, including both somatic and germ-line tissues.

Analysis of gene function through homology
The function of each gene was predicted based on the annotation categories assigned and their enrichment score in specific gene expression patterns observed in the beetle genome. This assignment allowed for the identification of over-represented gene functions present in each cluster using g: Profiler (P , 0.001) (Reimand et al. 2016) web server (Table S11 and S12).
Genes with functions such as cytoskeleton activity in cellular components (p-value = 7.8x10 28 ) and microtubule-based processes were enriched in Cluster-7 (Table S11). Genes associated with GO terms like transcription factor activity and binding activity as well as those that code for ubiquitin proteins, nucleus, nucleic acid metabolic process, DNA replication and microtubule motor activity proteins etc., were enriched in Cluster-10 (Table S11). Integral component of membrane was the most abundant cellular component category found in the ovary-expressed genes examined while in the biological process category; the majority of these genes were categorized as cell adhesion proteins (Table S11, Cluster-2). In the molecular function group, genes in Cluster-2 were enriched for ion binding and transporter proteins (Table S11).
Functional categories overrepresented in both germ-line and somatic tissues of male and female carcasses are presented in Table S12. In testis, enriched functional categories include microtubule-based processes; microtubule associated complex activities and movement of cell or subcellular component like biological process (Table S12). While the categories found highly enriched in genes present in both germ-cells include cell cycle (p-value = 2.1x10 210 ) and cell division (p-value = 2.9x10 24 ) (Table S12), genes with catalytic activity, metabolic process, hydrolase activity, binding activity, chitin and peptidases functions were found highly over-represented in soma-biased expressed genes (Table S12).

CONCLUSIONS
Important genetic and genomic resources, including a high quality genome sequence (Tribolium Genome Sequencing Consortium et al. 2008), have contributed to the ability to perform transcriptomic studies in Tribolium (Ninova et al. 2016;Preuss et al. 2012;Park et al. 2008;Altincicek et al. 2013), progressively making Tribolium a preferred model organism for molecular biology studies. In the present study, we have reported on the deep transcriptome of the Tribolium pre-zygotic embryo, zygotic embryo, testis, ovary, female adult soma, and male adult soma, and present the transcriptional profiles for three quarters (75%) of all predicted Tribolium genes in these tissues/stages (Tribolium Genome Sequencing Consortium et al. 2008).
Metazoan embryonic development is coordinated by both maternal RNAs deposited into the egg and the transcription of zygotic genes after fertilization. Depending on the insect species maternal transcripts persist a specific amount of time (Arbeitman et al. 2002). In Tribolium for example, we observed one set of maternal genes whose transcripts level substantially decreased after .6 hr following fertilization. and Enolase (C) was determined in testis tissue compared to ovary, male carcass, female carcass, embryo (0-5 hr) and embryo (6-11 hr) based on qPCR analysis. Fold change expression of four b-tubulins (Tc-b 1 , Tc-b 2 , Tc-b 3 and Tc-b 4 tubulin) in testis compared to ovary, male carcass, female carcass, embryo (0-5 hr) and embryo (6-11 hr) based on RNA seq analysis (D).
We observed the similarity and dissimilarity for maternally deposited homolog genes and early zygotic homologs in both Drosophila and Tribolium. The study of regulatory elements required for proper gene expression is crucial for understanding the molecular basis of evolution. Rather than the acquisition of new genes or modification of gene products, it is changes in gene regulatory elements that are often responsible for the evolution of new morphology (Carroll 2008). Wolff et al. (1998) studied both evolutionary conserved and diverged enhancers in Tribolium and Drosophila and reported that a major shift in the evolutionary transition from short (Tribolium) to long germ (Drosophila) embryogenesis is due to regulatory adaptation during the evolution of developmental processes (Wolff et al. 1998).
Gene expression levels are tightly regulated in ovary, testis, carcasses and early embryo and are likely controlled by cis elements that could be used to drive the expression of transgenic constructs for embryogenesis studies and genetic pest management strategies (Hammond et al. 2016;Buchman et al. 2018;Lai et al. 2018). Cis elements associated with genes identified in this study could be evaluated for functional activity using established reporter assays in Tribolium (Lai et al. 2018). Thus, the data we describe here can be mined for elements to construct synthetically engineered gene drive systems in Tribolium and to some extent other Coleopteran pest insects as well. Such gene drive systems have potential to replace, alter, or suppress wild populations of significant crop pests or other economical important insects (Buchman et al. 2018;Akbari et al. 2014). For example, recent strategies have featured the creation of male insects whose sperm express a transgene-based toxin that suppresses the insect population (Klein et al. 2012).

ACKNOWLEDGMENTS
This work was supported by the Defense Advanced Research Projects Agency project number HR0011-16-2-0036, as well as by Agrilife Research and the Department of Entomology at Texas A&M University. We thank David Pledger and Collin Valentin for assistance in insect rearing.

Figure 7
Highly expressed genes in the Tribolium ovary tissue. Venn diagrams representing increased gene expression greater than or equal to twofold in the Tribolium ovary in these comparisons (i.e., "Ovary vs Testis", "Ovary vs Male carcass" and "Ovary vs Female carcass").