Comparison of gene expression in the red imported fire ant (Solenopsis invicta) under different temperature conditions

The red imported fire ant (RIFA), Solenopsis invicta Buren is native to South America and is known as a global problematic invasive species. This study focused on the molecular response of RIFA by comparing gene expression profiles after exposing ants to low (10 °C) and high (40 °C) temperature stress and comparing them to untreated controls (30 °C). A total of 99,085 unigenes (the clustered non-redundant transcripts that are filtered from the longest assembled contigs) were obtained, of which 19,154 were annotated with gene descriptions, gene ontology terms, and metabolic pathways. 86 gene ontology (GO) functional sub-groups and 23 EggNOG terms resulted. Differentially expressed genes (DEGs) with log2FC ≥ 10 were screened and were compared at different temperatures. We found 203, 48, and 66 specific DEGs co-regulated at 10, 20, and 40 °C. Comparing transcriptome profiles for differential gene expression resulted in various DE genes, including cytochrome P450, NADH dehydrogenase subunit 1, cuticle protein and heat shock protein (HSP), which have previously been reported to be involved in cold and high temperature resistance. GO analysis revealed that antioxidant activity is up-regulated under high temperature stress. We verified the RNA-seq data by qPCR on 20 up- and down-regulated DEGs. These findings provide a basis for future understanding of the adaptation mechanisms of RIFA and the molecular mechanisms underlying the response to low and high temperatures.


Results
Sequencing, RNA-Seq assembly, and functional annotation. Quality filtering for Illumina raw data (Table S2) was completed to investigate the transcriptome responses to heat and cold stress in S. invicta. After transcriptome sequencing of four cDNA samples with Q30 > 94%, 44.53 GB of clean data passed the Illumina consistency filter (Table S3). All high-quality reads (Table S3) were pooled to perform the de novo transcriptome assembly. These contigs were further assembled into 107,264 transcripts with a mean length of 757.72 bp and a N50 of 1504 bp, and 99,085 unigenes with a mean length of 615.38 bp and a N50 of 1051 bp (Tables S4 and S5). The length distribution of unigenes was very similar to the transcript length distribution. This suggests a highquality assembly, which will serve as a sequence foundation for future research.
Annotation of predicted proteins. The assembled unigenes were validated and annotated using BLASTX against five public databases. Genes with a large blast hit to arthropods were detected after annotation. In total, 19,154 unigenes (19.33%) were discovered in at least one public database (UniProt). The NT database had the most matches (41,925 annotated unigenes, 42.31%), followed by the NR database (21,232,37.28%) (Fig. 1, Table 1). The majority of the unigenes were either unable to be annotated or had uninformative definitions (e.g., putative, unknown, hypothetical, or unnamed protein). According to BLASTX matches in the NR database, the unigene sequences were most similar to gene sequences from S. invicta (56.80%), and more than 70% showed similarity with ant genera (Solenopsis sp, Trachymyrmex sp, Acromyrmex sp, Atta sp, Camponotus sp, and Cyphomyrmex sp). ORFs with a duration of at least 100 amino acids were extracted. At least one ORF was found in 14.86% (14,721) of total expected unigenes (99,085), and 49.3% had a complete open reading frame ( Table 2).
GO and EggNOG analysis for global functional classification. The Gene Ontology (GO) database and the EggNOG database were used to identify the annotated unigenes. Since different functional annotations will exist for the same unigene, a total of 22,091 genes have been annotated to GO terms based on BLAST performance. Based on the three main GO groups, 'biological mechanism' , 'cell component' , and 'molecular function' the GO database yielded 86 GO functional sub-groups. The most frequent GO terms were 'metabolic process' , 'cellular process' , 'cell part' , 'catalytic activity' and 'binding' (Fig. 2A). The most common DEGs coregulated under cold stress according to GO categorization were 'cellular process' , 'metabolic process' , 'biological regulation' , 'developmental process' , 'cellular component' , 'organization or biogenesis' , and localization in biological processes in the biological process category. In the molecular function category, the coregulated DEGs were mostly assigned to 'cell par' , 'organelle' , 'membrane part' and 'Protein-containing complex' . For the cellular components category, only 'binding' and catalytic activity were significantly enriched (Table 3).

Differential gene expression under different temperatures.
We used the FPKM mapped reads method to measure the expression level of the unigenes. DEGs were found to be up-regulated under different temperature treatments as compared to untreated controls. With a criterion of p-value < 0.05 and |log 2 FC| ≥ 2, 4596, 2953, and 4068 unigenes were DEGs for T10, T20, and T40, respectively (Fig. 3A). A volcano plot with a criterion of p-value < 0.05 and log 2 FC ≥ 5 was constructed for each treatment temperature in comparison to the control temperature to identify the most likely temperature-reactive individual genes ( Fig. 3B-D). The largest number of expression changes in unigenes were observed in response to T10 (Fig. 3B). GO analysis revealed that most unigenes in the category of 'biological process' belong to the 'cellular process'  The X axis shows the ratio (%) of annotated unigenes for each database based on total unigenes. The Y axis shows the number of unigenes that are annotated based on each database.       (Table S6), and 'cellular part' and 'binding' included the largest number of unigenes from the 'cellular component' and 'molecular function' classifications, respectively (Table S6).
To explore the more specific and exclusive genes involved in cold and high temperature stress conditions, a Venn diagram was plotted for T10, T20, and T40 in comparison with T30 as a control group (p-value < 0.05). As shown in Fig. 4A, 203 (Table S9), 48 (Table S10), and 66 (Table S11) differentially expressed genes (DEGs) were identified when comparing the control (30 °C) and stressor temperatures of 10, 20, and 40 °C, respectively. There were 13 common DEGs that were consistently up-regulated in all three groups. The same 51 DEGs were up-regulated (FC ≥ 10) between T10 and T20, when compared with T30. There were 29 and 5 DEGs up-regulated (FC ≥ 10) between T10 and T40; T20 and T40, respectively, as compared to the T30 control. We believe two general groups of unigenes are related to temperature fluctuations. The first group includes the unigenes that are presented in all (T10, T20, and T40) or two (T10 and T40) treatment temperatures, and the second group includes the unigenes that are specifically expressed more than 10 times at one temperature (Fig. 4A). To understand the comparative distribution of unigenes in the first group, a heatmap was constructed (Fig. 4B,C). 'Venom carboxylase-6-like (p-value 3.29E−14)' , 'cGMP-dependent protein kinase (p-value 1.02E−09)' and 'growth hormone regulated TBC protein 1-A (p-value 4.36E−06)' showed high expression levels when the ants were incubated at 10 °C in comparison with 40 °C and 30 °C controls. Interestingly, 'histone' unigenes (histone H3-like centromeric protein (p-value 5.71E−06), histone H2A (p-value 1.21E−07), histone H4 (p-value 1.86E−06), and histone H2Blike (p-value 4.2E−05) showed high degrees of fold change at 40 °C in comparison with 10 °C (Fig. 4B). Among 13 matching unigenes, 'aromatic-l-amino-acid decarboxylase' and 'homeobox protein orthopedia-like' showed lower expression levels at T10 and T40 when compared to those at T20 (Fig. 4C, Table S8).
Specific unigenes (log 2 FC ≥ 10) for each treatment were the second group of unigenes that were investigated. KEGG and GO analysis were used to validate functional unigenes (Figs. 5 and 6, Tables S8, 9, and 10). KEGG pathway enrichment analysis revealed the primary DEG pathways ( Fig. 5A-C). Pathway enrichment was observed among all groups, although the T20 group showed less effects on pathway enrichment than T10 and T40 (Fig. 5B, Table S12 Table S12).
GO annotation was used to clarify the functions of DEGs that were significantly different (FC ≥ 10) between treatments (Fig. 6). Between T10 and T30 control, the most significant items were 'Biological Process: multicellular organismal process' , 'Cellular Component: organelle' and 'Molecular Function: catalytic activity and binding' (Fig. 6A). When comparing T20 and T30 controls, unigenes were assigned to 'Biological Process: developmental process' , 'Cellular component: membrane part' and 'Molecular Function: binding and molecular function regulator' (Fig. 6B). The comparison between T40 and T30 control revealed the following significant increase in unigene percentages: 'Biological Process: response to stimulus' , 'Cellular component: organelle' and 'Molecular Function: catalytic activity' (Fig. 6C).

Validation of gene expression profiles by q-PCR. qPCR and gel electrophoresis of twenty common
DEGs identified in the RNA sequence data were performed to confirm the accuracy and reproducibility of the Illumina RNASeq. The results of qPCR and Illumina FPKM ratio were plotted in Fig. 7. These data demonstrated that expression changes were in the same direction as the qPCR. The Illumina sequencing data was consistent with qPCR data, verifying the reliability and accuracy of the transcriptome analysis. This ensures the RNA-Seq results are considerably reliable for the identification of DEGs under temperature stress, and also the feasibility and sustainability of our further research into these or other DEGs from the transcriptome data.

Discussion
The red imported fire ant, Solenopsis invicta, is one of the most well-known ant invaders. The identification of geographical areas that are climatically vulnerable to S. invicta infestation should aid efforts to control the pest ant's spread 16 . S. invicta infestation threatens many parts of the world, including broad areas of Europe, Asia, Africa, Australia, and a number of island nations. Officials in charge of quarantine should be on the lookout for any unintentional introductions of this pest into vulnerable areas. As the size of the infestation increases, the cost of eradication rises sharply, and large infestations can be difficult to eradicate 16 . Researchers have long assumed that introduced S. invicta populations in the Southeast would not spread to higher latitudes and elevations because their pre-adaptation to the subtropical climate of their native range would render them unable    34 . With rising elevation, S. invicta's cold and heat tolerance thresholds decreased, suggesting a physiological capacity to withstand colder temperatures. The lipid content of S. invicta did not change with elevation, indicating that cold acclimation had no metabolic effects on the ants. These findings indicate that S. invicta is not as resistant to cold temperatures as previously believed, and that it will continue to adapt and spread to higher elevations and latitudes 35 . Comprehensive investigation of gene expression regulation under temperature stress is very important to understand the biochemical and physiological adaptation processes of invasive insect pests 22 . In this study, a comprehensive transcriptome analysis and characterization of the gene expression profiles of S. invicta under cold and high temperature stress were evaluated. Through the analysis of DEGs, transcriptome changes in S. invicta adult ants were revealed. Using RNA-seq techniques, four transcriptomes were de novo assembled from the adult stages of RIFA which were exposed to four different temperatures (10, 20, 30, and 40 °C), and 19,154 unigenes (19.33%) were successfully annotated from at least one public database (UniProt) ( Table 1). The results are in line with other transcriptome projects using Illumina technology [36][37][38] . More than 70% similarity with the ant genus was found and 56.80% of the unigene sequences were most similar to gene sequences from Solenopsis invicta. In this study, DEGs from adult RIFA subjected to various treatment temperatures (10, 20, and 40 °C) were compared to a 30 °C control group in this analysis. The majority of DEGs were observed at T10, followed by T40, both of which expressed a greater DEG distribution than T20 (Fig. 3A). As mentioned earlier, this is consistent with proteomics data from Locusta migratoria under high and low temperature stress 39 . To identify specific genes associated with response to temperature, the number of unigenes with log 2 FC ≥ 10 was clarified by the Venn diagram, and KEGG analysis www.nature.com/scientificreports/ was conducted to determine the probability of function in pathway enrichment. KEGG analysis revealed that of 203 specific cold-regulated DEGs (associated with T10), 41 DEGs were enriched in the KEGG pathways; 'Metabolic pathway' , 'Carbon metabolism' , 'Citrate cycle (TCA)' , 'RNA transport' , and 'Lysosome' . Interestingly, in T20 and T40, the 'Metabolic pathway' included more DEGs than other pathways (Fig. 5). 'Purine metabolism' , 'Spliceosome' , 'Lysosome' , 'RNA degradation' , 'Glycolysis/Gluconeogenesis' , 'Pyruvate metabolism' , 'Phagosome' , 'Sphingolipid metabolism' , 'RNA transport' , 'Glycerolipid metabolism' , 'Carbon metabolism' , 'ECM-receptor interaction' are pathways that demonstrate similar results as those of investigations on transcriptome responses to cold stress in the carpenter moth, Eogystia hippophaecolus 36 , and the a chrysomelid beetle, Galeruca daurica 22 .
The transcriptome analysis revealed that the expression of 'Glycolysis' and 'TCA cycle pathways' are up-regulated in a similar manner to the braconid wasp, Aphidius colemani, when exposed to low temperatures 40 . When RIFA was exposed to the highest temperature (40 °C), 'Tyrosine metabolism' , 'Phenylalanine metabolism' , 'Cysteine and Methionine metabolism' , 'Spliceosome' , 'Protein processing in endoplasmic reticulum' , and 'Metabolic pathway' , were found to be enriched pathways that are similarly enriched in three species' of rice plant hopper when exposed to 37 °C 28 . 'Fatty acid synthase' and 'Fatty acid metabolism' , are two of the main pathways of RIFA when exposed to high temperatures. The expression of fatty acids as hydrophobic agents allows insects to avoid water loss in warmer regions of the globe 41 . Gomphocerus sibiricus is known to increase its levels of oleic acid, linoleic acid, and glycerin in response to high temperatures, a process that may reduce mortality due to excessive evaporation of body moisture 42 . In our study on RIFA, ' Amino acid metabolism' was clearly up-regulated during high-temperature stress. It is suggested that amino acid metabolism provides heat resistance in RIFA similar to that of results that have been reported for Locusta migratoria 39 when exposed at 40 °C. Due to the synthesis of immune proteins and defense enzymes, insects seek out and consume numerous free amino acids when coping with stress conditions such as high temperature, low temperature, and fungal invasion 43 . The synthesis and metabolism of amino acids are necessary to produce a significant number of amino acids, which make available the raw materials necessary for the synthesis of heat-resistant proteins 39 .
In this study, two cuticular protein unigenes were identified from 203 co-regulated DEGs under cold temperature stress (Table S9). Cuticular protein gene expression has been observed in studies of other insects such as beetles, moths, planthoppers, and stick insects when exposed to cold temperature stress 22,28,36,44 . Although the physiological role of cuticular proteins in insect cold hardiness has not yet been identified, it seems insect cuticle may play an important role in insects when coping with low temperatures 22,25,28,45,46 .
According to GO analysis (Fig. 6A), ' Antioxidant activity' was enriched at low temperatures suggest that it might contribute to RIFA ability to resist oxidative stress damage at low temperatures 22 , or their potential for cell preservation via antioxidant defense when challenged by environmental complexity 47 . In addition, one cytochrome P450 was identified that is up-regulated exclusively under low temperatures (Table S9). Meanwhile, NADH dehydrogenase subunit 1 was up-regulated at all treated temperatures (Fig. 5C, Table S8). These two proteins are the main enzymes in the antioxidant activity pathway 36 . In a comparative analysis of the transcriptional responses to low and high temperatures in three rice planthopper species, some cytochrome P450 genes were up-regulated under both low and high temperatures, which suggests cold and heat stress increase oxidative stress in the insect body 28 .
Heat shock proteins (HSPs) are another important protein that insects use as critical physiological products when under abiotic stress conditions 39 . From earlier studies, it was believed that Hsp is associated with biological cold and heat resistance 48,49 . HSPs are molecular chaperones, that play important physiological roles, including correct folding of proteins, prevention of protein denaturation, and degradation of misfolded or condensed proteins, and maintenance of correct protein conformation 50,51 . In this study, we identified one specific Hsp70 that up-regulated about a 12-fold change when RIFA was exposed at 10 °C (Table S9). Two pathways, including 'Protein export' and 'Protein processing in endoplasmic reticulum' , were enriched under Hsp70 gene expression at low temperatures (Fig. 5A, Table S12). Interestingly, heat shock protein 83 was found only in the 40 °C treatment group that was up-regulated about 11 times (Table S11). Many studies have confirmed that the expression of Hsp genes can be up-regulated by cold and heat stimulus 50,52 . To assist the resistance to temperature stress, the Hsp60 gene expression in Stegobium paniceum significantly increases under high-and low-temperature stress 53 . Three Hsp90 and four Hsp70 were up-regulated by cold stress and were differentially expressed in the desert beetle, Microdera punctipennis 27 . The differences in Hsp, insect species, sex of organism, and intensity of temperature are important factors related to Hsp expression level in insects 22,54 .
In conclusion, we compared the transcriptomes of S. invicta under high-and low-temperature stresses using RNA-Seq technology based on high-throughput sequencing. Comparative transcriptome analysis identified many genes, and a large number of changes were discovered in metabolic pathways through GO and KEGG enrichment analysis. Our data will facilitate further molecular investigations and genomic research. Many novel relationships between high-and low-temperature and significantly up-regulated genes were identified in this study (Tables S7-11). These newly found genes may be important for RIFA overwintering and adaptation potential in new environments as well as quarantine areas.

Materials and methods
Insect rearing, exposure temperatures and sample preparation. Solenopsis p  s  e  p  a  rt  c  e  ll  p  a  rt  s  y  n  a  p  s  e  n  u  c  le  o  id   s  u  p  ra  m  o  le  c  u  la  r  c  o  m  p  le  x   U  n  c  la  s  s  if  ie  d  2   c  a  ta  ly  ti  c  a  c  ti  v  it  y   s  ig  n  a  l  tr  a  n  s  d  u  c  e  r  a  c  ti  v  it  y   s  tr  u  c  tu  ra  l  m  o  le  c  u  le  a  c  ti  v  it  y   tr  a  n  s  p  o  rt  e  r  a  c  ti  v  it  y  b  in  d  in  g   a  n  ti  o  x  id  a  n  t  a  c  ti  v  it  y   p  ro  te  in  ta  g   c  a  rg  o  re  c  e  p  to  r  a  c  ti  v  it  y   tr  a  n  s  la  ti  o  n  re  g  u  la  to  r  a  c  ti  v  it  y   m  o  le  c  u  la  r  fu  n  c  ti  o  n  re  g  u  la  to  r   m  o  le  c  u  la  r  c  a  rr  ie  r  a  c  ti  v p  s  e  p  a  rt  c  e  ll  p  a  rt  s  y  n  a  p  s  e   s  u  p  ra  m  o  le  c  u  la  r  c  o  m  p  le  x   U  n  c  la  s  s  if  ie  d  2   o  b  s  o  le  te  R  N  A  p  o  ly  m  e  ra  s  e  II  I  ty  p  e  3  p  ro  m  o  te  r  T  F  II  IB  -t  y  p  e  T  F  a  c  ti  v  it  y   c  a  ta  ly  ti  c  a  c  ti  v  it  y   s  ig  n  a  l  tr  a  n  s  d  u  c  e  r  a  c  ti  v  it  y   s  tr  u  c  tu  ra  l  m  o  le  c  u  le  a  c  ti  v  it  y   tr  a  n  s  p  o  rt  e  r  a  c  ti  v  it  y  b  in  d  in  g   a  n  ti  o  x  id  a  n  t  a  c  ti  v  it  y   p  ro  te  in  ta  g   c  a  rg  o  re  c  e  p  to  r  a  c  ti  v  it  y   m  o  le  c  u  la  r  fu  n  c  ti  o  n  re  g  u  la  to  r   m  o  le  c  u  la  r  c  a  rr  ie  r  a  c   www.nature.com/scientificreports/ trays so that any escaping ants would be trapped inside the holding trays. Ants were fed a 20% sucrose solution, mealworms (Tenebrio molitor larvae), and an artificial diet described by Dussutour and Simpson 55 . Water was provided ad libitum. Colonies contained dealated mated queens, alate queens, males, brood (eggs, larvae, and pupae) and workers. To perform transcriptomic analysis, medium sized workers were incubated at 10, 20, and 40 °C for 24 h. Ants incubated at 30 °C were considered as the control group. A temperature recorder (Lutron, BTM-4208SD) was used to check the output temperature, which was recorded every 2 s for 24 h. According to the collected data, the incubator's variance is 0.5°. Each treatment was replicated three times. After the temperature treatment, ten ants from each group were immediately frozen in liquid nitrogen and stored at − 80 °C for subsequent experiments.
De novo assembly. Illumina short reads were quality-filtered and adapter-trimmed using Trimmomatic v0.38 (http:// www. usade llab. org/ cms/? page= trimm omatic). FastQC v0.11.7 (http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/) was used to check data quality before and after trimming. After the removal of lowquality reads, an Illumina-based de novo transcriptome assembly was performed using Trinity version trinity rnaseq r20140717, bowtie 1.1.2 58 . Trimmed reads for every sample were merged into one file to construct a combined reference. The de novo assembly of merged data was carried out using Trinity with default parameters and assembled into transcript contigs 59 . The total number of genes, transcripts, GC content, max/min/median/ average contig length, and total assembled bases were summarized. Trinity groups transcripts into clusters based on shared sequence content. For assembled genes, the longest contigs of the assembled contigs are filtered and clustered into non-redundant transcripts using CD-HIT version 4.6 (http:// weizh ongli-lab. org/ cd-hit) 60 . These transcripts were defined as 'unigenes' which are used for predicting ORFs (Open Reading Frames), annotating against several known sequence databases, and analyzing differentially expressed genes (DEGs Differential gene expression analysis. A quality check was conducted for all samples, so that if more than one read count value was zero, it was not included in the analysis. Gene expression levels were measured in the RNA-Seq analysis as fragments per kilobase of transcript per million mapped reads (FPKM) 69 . Multiple testing was corrected for in all statistical tests using the Benjamini-Hochberg false discovery rate with the following parameter values: FDR < 0.01 36 . In order to reduce systematic bias, we estimated the size factors from the count data and applied Relative Log Expression (RLE) normalization with the DESeq 2 R library. Using each sample's normalized value, the high expression similarities were grouped together by Hierarchical Clustering Analysis and graphically shown in a 2D plot to show the variability of the total data using Multidimensional Scaling Analysis. Significant unigene results were analyzed as Up and Down-regulated count by log 2 FC ≥ 5, ≤ − www.nature.com/scientificreports/ 5 and ≥ 10, distribution of expression levels between the two groups was plotted as Volcano plot (https:// huyge ns. scien ce. uva. nl/ Volca NoseR) 70 and simple bar plots. The DEGs were then used for GO and KEGG enrichment analysis using the edgeR exact test. The software topGO was used to carry out GO enrichment analysis. All DEGs were aligned to terms in the KEGG database and searched for significantly enriched KEGG terms. Heat maps were generated using the online tool Heatmapper (http:// www. heatm apper. ca/ expre ssion/) 71 .