In silico-driven analysis of the Glossina morsitans morsitans antennae transcriptome in response to repellent or attractant compounds

Sample collection, RNA extraction and sequencing

The male G. m. morsitans flies used in this study were from a colony reared at Yale University insectary as described by Bateta et al. (2017). Feeding, exposure to the chemicals and RNA extraction was done as described (Kabaka et al., 2020). A pair of antennae were carefully hand-dissected from fifty flies in each treatment (attractant, repellent, or control) and replicate (1, 2 or 3) as described (Menuz et al., 2014). PCR amplicons were generated using tsetse fly specific beta-tubulin gene primers to confirm removal of the gDNA from total RNA (Bateta et al., 2017). Sequencing was done on Illumina HiSeq 2500 at Yale University Center of Genome Analysis (YCGA), New Haven, CT, USA. The raw transcriptome sequences have been deposited at the Sequence Read Archive (SRA) under study accession number PRJNA343267.

Data preprocessing and differential gene expression analysis

Quality of the FastQ files was checked using FastQC v0.11.5 (Babraham Bioinformatics, 2013) software and the output cleaned using Trimmomatic software v0.38 (Bolger, Lohse & Usadel, 2014). Any contaminating rRNA reads in the data were removed using SortMeRna v2.0 (Kopylova, Noé & Touzet, 2012). The clean paired reads from different treatments and replicates were then separately mapped onto G. m. morsitans transcripts gene-set version 1.9 from Vectorbase or genome version 1.0 (Giraldo-Calderón et al., 2015) using STAR software v2.7.3a (Dobin et al., 2013). VectorBase provides reliable datasets that are community reviewed and regularly updated. Chemosensory proteins in these datasets have been annotated (Obiero et al., 2014; Macharia et al., 2016; Liu et al., 2012). Transcript mapping provided information on transcript specific mapping abundance (Mortazavi et al., 2008). Mapping reads to the genome was an additional quality control procedure that would get rid non-G. m. morsitans contaminants. The BAM files contained information on sequences aligned onto the transcripts, sorted by respective coordinates to facilitate downstream analyses. The number of reads (counts) aligned onto each transcript in the respective BAM files were quantified using Salmon software v1.2.1 (Patro et al., 2017). This analysis provided data on relative abundance of read from different treatments and replicates that mapped onto G. m. morsitans transcripts.

Differential expression analysis was done on the gene count matrix from Salmon using the DESeq2 R-package version 1.28.0 (Love, Huber & Anders, 2014). Default parameters for count data normalization as recommended (Conesa et al., 2016) were used to allow for control of log₂ fold change shrinkage, custom p-value and fold change cut-offs. Genes were considered differentially expressed and retained for further analysis if the test statistics p-value (adjusted for false detection rate) (FDR) was less than 0.05 according to the method from Benjamini & Hochberg (1995). Since the antennae is functionally specialized for olfaction, and potentially enriched with associated canonical chemosensory gene transcripts (Kabaka et al., 2020), we separately probed for expression profiles of these transcripts, and isolated those with at least two-fold change in difference between the attractant or repellent and control treatments. We generated heatmaps for visualizing the relationships between the different treatments using the package Pheatmap v1.0.12 (Kolde, 2012) in R software (R Core Team, 2017).

Co-expression analysis using WGCNA

WGCNA package is designed for clustering genes based on their expression profiles and therefore we used the gene lists generated from Salmon for weighted co-express analysis. Normalization was done by filtering out genes with counts less than 10 in more than 90% of samples since they are not informative and tend to introduce noise (Farhadian et al., 2021). Co-expression networks were generated using WGCNA using the Bioconductor R package v3.5.1 (Langfelder & Horvath, 2008). This algorithm calculates a similarity co-expression matrix using correlation for all genes defined as: $s_{ij}=\left|\mathrm{c}\mathrm{o}\mathrm{r}\left(x_{i_{}},x_j\right)\right|$. An adjacency matrix was calculated by raising the co-expression similarity to a soft thresholding power beta $\beta$ defined as $a_{ij}=p_{ower(S_{\mathrm{i}\mathrm{j}}},\beta )={\left|s_{ij}\right|}^{\beta }$, where $a_{ij}$ represents the resulting adjacency which is a measure of the connection strengths. A $\beta$ of 12 was selected for construction of a gene co-expression network based on estimated scale-free topology as described (Langfelder & Horvath, 2008). A topological overlap matrix (TOM) was computed, converted into dissimilarity matrix and then hierarchical clustering used to generate a tree (dendrogram). Co-expression modules were detected using Dynamic Tree Cut (DTC) algorithm based on an edge height cut-off of 0.25, and a minimum module size of 30 genes. Relationship between different co-expression was explored using FlashClust function (Langfelder & Horvath, 2012) and a heatmap was used to visualize the correlation between the modules.

Network analysis

Network visualization and centrality measure analysis were done using Cytoscape v3.7.2 (Shannon et al., 2003). The degree, betweenness centrality and clustering coefficient of the network were analyzed using network analyzer (a Cytoscape plugin) as described by Assenov et al. (2008). Degree $d\left(v\right)$ of a vertex v, in a network $G=\left(V,E\right)$, counts the number of edges in E incident upon v. Given G, define f (d) to be the fraction of vertexes $v\in V$ with degree $d\left(v\right)=d$. Genes having large degrees are referred to as hubs genes, indicating that they hold multiple genes/proteins together and have the highest potential to regulate the node $v$. Betweenness is a measure of the number of shortest paths that are connecting any two nodes (j, k). On the other hand, closeness is a measure of the ability of a node to interact with all other nodes, including the indirectly connected nodes. It is defined as:

$$\sum _J^{j_1={\displaystyle \frac{1}{N}}}h\left(\mathrm{i},i\right)$$

with $h\left(\mathrm{i},i\right)$ being the shortest distance between gene $\mathrm{i}$ and $\mathrm{j}$. Finally, clustering coefficient is defined as the edge density of the neighborhood of node $i$ which is calculated as:

$$u_{i={\displaystyle \frac{m_i\left(j,k\right)}{m_i}}}$$

where $m_i\left(j,k\right)$ is the number of edges connecting nodes (j, k) neighboring node i and m_i is the total number of visible edges of all the neighboring nodes of i that are fully connected.

Discretization and association rule mining

Discretization is a data pre-processing step used in machine learning in order to transform continuous or numerical attributes into discrete ones (Kotsiantis & Kanellopoulos, 2006; Hacibeyoglu & Ibrahim, 2018). In this study, we used the equal frequency discretization method (Dougherty, Kohavi & Sahami, 1995) implemented in GEDPROTOOLS (Gallo et al., 2016) to transform the data from continuous to discrete values. Equal frequency discretization reduces the effect of outliers and collects similar values in the same interval (Hacibeyoglu & Ibrahim, 2018). Genes that were co-expressed in the network were retrieved for discretization using the steps outlined below:

• Input: the continuous values of attribute and number if intervals A = {a₁, a₂,…, a_n−1, a_n}

  and number of intervals k, where k > 0.

  Step 1: Sort all values of in ascending order,

  Step 2: Divide A by k intervals,

  Step 3: Create bins according to number of elements in each interval,

  Step 4: Determine boundaries of each interval by calculating the average value of the

  Maximum value of the current bin and the minimum value of the next bin,

  Step 5: The continuous values of A are transformed into discrete ones by determining

  the interval that they belong to,

  Output: A with discrete values

Two bins were created in step 3, with the final output being discretized measurements whereby a value of zero represented a gene that was under-expressed, while a value of one represented a gene that was overexpressed. In this context, the discretized data was used to identify frequent itemsets using the Apriori algorithm (Agrawal & Srikant, 1994) implemented in the R package arules ver. v1.6-4 (Hahsler et al., 2014). In association rule mining, a rule is typically described by three measures: support, confidence, and lift. These three represent the significance and interest of a rule. Support of a rule $X\Rightarrow Y$ is equal to the support of the itemset $X\cup Y$ and is defined as the probability of finding all the genes in sets X. Support of an itemset X is calculated as:

$$Suppor{t}_D\left(x\right)={\displaystyle \frac{\left|\left\{T\in Dx\subseteq T\right\}\right|}{\left|D\right|}}$$

The confidence of rule $X\Rightarrow Y$ is the probability of finding all the differentially expressed genes in set Y as compared with the differentially expressed genes in set X. The confidence is calculated as:

$$Confidenc{e}_D\left(x\Rightarrow y\right)={\displaystyle \frac{Sup{p}_D\left(X\cup Y\right)}{sup{p}_D\left(X\right)}}$$

Lift measures the strength of the rule and varies in the interval [0, ∞]. Lift is defined as:

$$lift\left(X\to Y\right)={\displaystyle \frac{supp\left(X\cup Y\right)}{supp\left(X\right)\ast supp\left(Y\right)}}$$

We used a support value of 0.5 and a confidence value of 0.99. Support values greater than 0.5 gave zero rules while support values of less than 0.5 and confidence of less than 0.99 resulted in too many rules. We therefore filtered and retained only the rules that had a lift value ≥ 2.

Gene set enrichment analysis

Gene Set Enrichment Analysis (GSEA) was done using WEB-based Gene SeT AnaLysis Toolkit (WebGestalt) as described by Wang et al. (2017). Using Drosophila melanogaster homologs as proxy to assess gene enrichment since the database of WebGestalt does not have G. morsitans, homologs. FDR corrected and p-value ranked D. melanogaster gene homologs of the differentially expressed G. m. morsitans genes were selected as input for the analysis using default parameters (5–2,000 Entrez Gene IDs, FDR < 0.05, 1,000 permutations and 20 categories with the outputted leading-edge genes) as described by Dębski et al. (2016). We separated significantly enriched non-redundant biological processes, cellular components, and molecular function Gene Ontology (GO) terms as identified by GSEA (Ashburner et al., 2000; Stark et al., 2006; The Gene Ontology Consortium, 2017).

Differential expression in response to repellants (δ-nonalactone) or attractant (Ɛ-nonalactone)

Raw reads per sample ranged between 23 and 73 million. Ribosomal RNA contamination levels were below 5% in all the samples indicating a near-perfect ribosomal RNA depletion before sequencing was done. The number of clean reads in the samples after preprocessing and quality filtering ranged between 11.1 M and 38.3 M. Reads that mapped to the genome were between 78% and 92%, while the remainder mapped to multiple locations, were not unique or had no features. During differential gene expression analysis, 2,097 low-count genes were filtered out leaving a total 10,921. Genes with a False discovery rate (FDR) < 0.05 were considered as being of biological significance (Von Der Weid et al., 2015). A set of 308 genes were identified as differentially expressed after exposure to either repellent (δ-nonalactone) or attractant (Ɛ-nonalactone) as compared to the control (no treatment). Figure 1A shows the top 50 significantly expressed chemosensory genes (FDR < 0.05) in response to Ɛ-nonalactone. Nine genes showed upregulation in flies exposed to an attractant (Ɛ-nonalactone). These are the Gustatory receptors (GrIA and Gr28b), ionotrophic receptors (Ir41a and Ir75a), odorant binding proteins (Obp99b, Obp99d, Obp59a and Obp28a) and the odorant receptor (Or67d). In addition, several non-chemosensory genes with no assigned function in the NCBI database were co-expressed with the chemosensory genes. In contrast, exposure to repellent (δ-nonalactone) did not show any significant change between the treatment and control samples. We identified 24 non-chemosensory that showed significant upregulation in response to Ɛ-nonalactone as shown in Fig. 1B. Six of these genes have no assigned function biological function and were also upregulated in the chemosensory geneset.

Co-expression network analysis and hub genes identification

A scale-free topology weighted gene network was constructed using WGCNA based on a soft thresholding power (β). From candidate powers of between 1–20, β = 12 returned a scale-free topology fit index of –0.1. An adjacency matrix based on the criterion of approximate scale-free topology is shown in Figs. 2A and 2B. Using the dynamic tree cutting algorithm, all the 11,808 genes were grouped into 12 modules, which ranged in size from 42 to 9,325 genes per modules (Figs. 2C and 2D).

We identified 12 modules (Fig. 3) with the major modules being turquoise (n = 9,325 genes), blue (n = 1,040 genes), and brown (n = 377 genes).

Correlation between modules

We performed cluster analysis to identify whether chemosensory genes were evenly distributed with the 12 co-expressed modules. Using an edge height cut-off of 0.25, all modules were below the 75% similarity hence no modules were merged. Interestingly, most of the chemosensory genes (n = 83) were in the turquoise module which also had 77 non-chemosensory genes as well as 46 genes that have no assigned function yet. In the other modules, genes Or67d was identified in the blue module, Obp19b & Obp83g in the brown module, and Obp73a in the yellow module. Therefore, the turquoise module was interesting from the chemosensation point of view. We proceeded to probe if there exist relationships between genes in the turquoise, blue, brown and yellow modules by visualizing only the functionally annotated genes in a network (Figs. 4A and 4B). Chemosensory genes are depicted in green, non-chemosensory in yellow and those with unknown function in grey. We filtered out genes with degree value of less than 5 to reduce the size of the graph and got a final network with 51 nodes and 148 edges (Fig. 4B). Filtering changed the network density from 0.05 to 0.116 while overall clustering coefficient increased from 0.475 to 0.587 (Figs. 4C and 4D).

The degree, average shortest path length and clustering coefficient for the top genes with a node greater than 8 are shown in Table 1. Genes CAH3, Ahcy, Ir64a, Or67c, Ir8a and Or67a can be regarded as the top hub genes since they had degree values above 10. Fourteen of the top 20 genes are associated with chemosensation, which is an important biological function in insects. Clustering coefficient of hub gene nodes ranged between 0.29 and 0.68 and this is an indication that some parts of the network were more intricately connected than others.

Association rule mining

Discretized and transformed gene expression data is presented in a transaction format where the samples represent transaction IDs and genes represent items. A minimum support of 0.5 and a confidence of 0.99 was used to generate 801 rules, which we further filtered using a lift of $\ge 2$ to empirically produce the best results. Lift values lower than 2 or support values less than 0.5 generated too many rules whereas when we used support values greater than 0.5, no rules were generated. Genes with no assigned biological function were of interest since we wanted to find out if association rule mining could help predict their function. We therefore arrowed down the rules that implied an association between known genes and those with no known function. Twenty-two representative rules are shown in Table 2.

Gene set enrichment analysis (GSEA)

GSEA of the antennae transcripts revealed several enriched processes involved in detection of stimulus. GoSlim GO analysis component of the GSEA assigned most of the significantly expressed genes to various biological, cellular components and molecular function as shown in Fig. 5.

The top Enriched biological processes included biological regulation, response to stimulus, metabolic processes, and cell communication. Enriched cellular components were mainly associated with the membrane. Most enriched molecular functions were associated with ion binding, molecular transducer activity and protein binding. As we expected, overrepresentation analysis as well as the network expansion analysis showed that most of the significantly expressed genes were associated with sensory responses.