Comparison of Gene Expression in the Red Imported Fire Ant, Solenopsis Invicta, in Different Temperature Conditions

The red imported re 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 proles after exposing ants to low (10 ℃ ) and high (40 ℃ ) temperature stress and comparing to untreated controls (30 ℃ ). A total of 99,085 unigenes 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 log 2 FC ≥ 10 were screened and were compared at different temperatures. We found 203, 48, and 66 specic DEGs co-regulated at 10, 20, and 40 ℃ . Comparing transcriptome proles for differential gene expression resulted in various DE proteins and 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 up-regulated under high temperature stress. We veried the RNA-seq data by qPCR on 20 up and down-regulated DEGs. These nding provide a basis for the future understanding of adaptation mechanism of RIFA and molecular mechanism underlying the response to low and high temperatures. dT to the All quantitative (qPCRs) using a Real-Time The forward and µL nuclease free water. RT-qPCR cycling with 95°C treatment followed by at for at for and for The expression level of as a gene used to normalize target gene expression levels 45 under different treatments. PCR products assessed by melting curve analysis. Quantitative analysis performed using comparative −ΔΔCT


Introduction
The red imported re ant (RIFA), Solenopsis invicta Buren (Hymenoptera: Formicidae), is a global invasive and aggressive species native to South America that is considered to be one of the 100 most impactful invasive pests around the world 1,2 . RIFA was rst reported from the US in the 1930's and since has spread to other temperate areas around the globe 2011 1 . RIFA populations are now established in the United States, Mexico, Australia, New Zealand, China, Malaysia, Japan, Singapore, and the West Indies 3,4

. In
Korea, RIFA is listed as a quarantine pest. RIFA is well-documented to drive negative impacts human health, public safety, ecosystems, agriculture and native biodiversity in their invasive range 5 . RIFA invasions potentially threaten 41 species on the China National List of Protected Wildlife, including 22 birds, one amphibian and 18 reptiles and also it has been predicted that they will likely create resource limitations in arthropod communities 6 . As a result of the deleterious effects of RIFA invasions mentioned above, RIFA is of great concern to the Animal and Plant Quarantine Agency of Korea.
Temperature is one of the critical abiotic factors that determine the distribution and life history of insects 7,8 . Population dynamics and geographical distribution of insects are affected by temperatures trough interfering with their metabolic processes such as alimentation, digestion, detoxi cation, mating, and development [9][10][11][12][13] . Extreme temperatures are potential hazards for the stability of insect populations and can be detrimental to their development. In warm regions, RIFA demonstrates ecological adaptability to extreme high temperatures, but its geographical distribution is directly impacted by cold temperatures 14 . It has been well documented that temperature indices provide a useful predictive tool for predicting the potential distribution of RIFA in newly invaded systems 15 . RIFA is of great concern in China, where studies have been conducted to determine its tolerance to extreme temperatures in order to predict its potential range expansion 14 . Recently, scientists using 'omic' technologies have determined which pathways are important for allowing a species of beetle to cope in temperature stress 16 . Transcriptomics and the fast development of novel high-throughput sequencing technologies, such as RNA-Seq, has provided an opportunity to investigate signaling-associated genes and triggered putative function(s) and pathway(s) at low and high temperature stress conditions, in insects 17,18 .
RNA-seq technologies in a New Zealand alpine stick insect demonstrated upregulation of cuticle genes following cuticle modi cation in response to low temperature were the results of rst used 19 . Since 2014, transcriptome analysis using RNA-seq has been used to investigate gene expression changes when coping with thermal stress in several species of insects (Drosophila virilis 20 , Cryptolaemus montrouzieri 18 , Microdera punctipennis 21 , Nilaparvata lugens, Sogatella furcifera, Laodelphax striatellus 22 , Galeruca daurica 16 , and Monochamus alternatus 7 ). The ndings of such studies demonstrate that cold stress can change the expression levels of hundreds of genes associated with transcription, metabolism, and cuticular organization, especially enzyme-related genes responsible for the upregulation of encoding cytochrome P450s (P450), antioxidative enzymes, and aldehyde dehydrogenase 18,23,24 .
In this study, we used RNA-Seq and de novo transcriptome assemblies to generate transcriptomes and examine the changes in the regulation of transcription associated with cold and heat treatment in S. invicta. Detailed differential expression analysis revealed a number of candidate genes that are potentially related to the cold and heat tolerance of RIFA. We performed qRT-PCR to validate the RNA-seq data. We aimed to develop a basis for the adaptive mechanism and a rich resource for the discovery and identi cation of novel genes involved in the cold and heat stress response in S. invicta.

Sequencing, RNA-Seq Assembly, and Functional Annotation
To investigate the transcriptome responses to heat and cold stress in S. invicta, quality ltering for Illumina raw data (Table S2) was completed. In total, 44.53 Gb of clean data passed the Illumina quality lter after transcriptome sequencing of four cDNA samples with Q30 > 94% (Table S3). To perform the de novo transcriptome assembly, all high-quality reads (Table S3) were pooled. Using paired-end joining and clustering according to the similarity of contigs, these contigs were further assembled into 107,264 transcripts with a mean length of 757.72 bp and an N50 of 1,504 bp, and 99,085 unigenes with a mean length of 615.38 bp and an N50 of 1,051 bp (Table S4 and S5). The length distribution of unigenes closely followed the length distribution of transcripts. This indicates a high-quality assembly, providing a sequence basis for future studies.

Annotation of predicted proteins
For validating and annotating, the assembled unigenes were searched against ve public databases (NR, NT, UniProt, Pfam, GO, EggNOG, and KEGG) using BLASTX with a cut off E-value of 10 − 5 . After annotation, genes with a signi cant blast hit to arthropods were identi ed. In total, 19,154 (19.33%) unigenes were found in at least one public database (UniProt). The NT database (41,925 annotated unigenes, 42.31%) had the most matches, followed by the NR database (21,232, 37.28%) (Fig. 1, Table 1).
Overall, most of the unigenes either could not be annotated or their descriptions were uninformative (e.g., putative, unknown, hypothetical, or unnamed protein). Overall, 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), as observed via BLASTX matches in NR database.  GO functional sub-groups could be obtained according to the three main GO groups 'biological process', 'cell component', and 'molecular function' in the GO database. 8,646 unigenes belonged to the biological process group, 6,676 unigenes fell in the cellular component group, and 6,769 unigenes were categorized in the molecular function group ( Fig. 2A). The most frequent GO terms were 'metabolic process (9,546 unigenes)', 'cellular process (12,843 unigenes)', 'cell part (14,118 unigenes)', 'catalytic activity (8,895 unigenes)' and 'binding (9,186 unigenes)'.
We used the EggNOG database to reveal functional and biological classi cation. In total, 31,092 unigenes were assigned to 23 EggNOG terms (Fig. 2B) that belonged to three functional classes including 'information storage and processing', 'cellular processes and signaling', and 'metabolism'. The largest number of unigenes were classi ed as 'translation, ribosomal structure, and biogenesis (
To identify the most probable temperature-reactive speci c genes, a volcano plot with a criterion of pvalue < 0.05 and log 2 FC ≥ 5 was constructed for each treatment temperature, in comparison as the control temperature ( 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 (12,843)' sub-group (Table S6), and 'cellular part (14,118)' and 'binding (9,186)' included the largest number of unigenes from the 'cellular component' and 'molecular function' classi cations, respectively (Table S6).
To explore the more speci c and exclusive genes involved in cold and high temperatures stress conditions, a Venn diagram was plotted for T10, T20, and T40 in comparison with the T30 as a control group (p-value < 0.05). As shown in Fig. 4A, 203, 48, and 66 differentially expressed genes (DEGs) were identi ed when comparing the control (30°C) and stressor temperature 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 temperatures uctuations. The rst group includes the unigenes that are presented in all (T10, T20, and T40 (13)) or two (T10 and T40 (42)) treatment temperatures, and the second group are the unigenes that speci cally expressed more than 10 times at one temperature (T10 (203), T20 (48), and T40 (66)) ( Fig. 4A). To understand the comparative distribution of unigenes in the rst group, a heatmap was constructed ( Fig. 4B and C). 'Venom carboxylase-6-like', 'cGMP-dependent protein kinase' and 'growth hormone regulated TBC protein 1-A' showed high expression levels when the ants were incubated at 10°C in comparison with 40°C and to 30°C controls. Interestingly, 'histone' unigenes (histone H3-like centromeric protein, histone H2A, histone H4, and histone H2B-like) showed high degrees of fold change at 40°C in comparison with 10°C (Fig. 4B).

Discussion
Comprehensive investigation of gene expression regulation under temperature stress is very important to understand the biochemical and physiological adaptation processes in invasive insect pests 16 . In this study, a comprehensive transcriptome analysis and characterization of the gene expression pro les 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 was de novo assembled from the adult stages of RIFA which 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 databases (UniProt) ( Table 1). The results are in line with other transcriptome projects using Illumina technology [25][26][27] . 56.80% of the unigene sequences were most similar to gene sequences from Solenopsis invicta and more than 70% similarity with ant genus were observed. In this study, DEGs from adult RIFA which were exposed to different treatment temperatures (10, 20, and 40°C) were compared to that of a 30°C control group. The majority of DEGs were observed at T10, follow 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 28 . To identify speci c genes associated with response to temperature, the number of unigenes with log 2 FC ≥ 10 was clari ed by Venn diagram, and KEGG analysis was conducted to determine the probability of function in pathway enrichment. KEGG analysis revealed that from 203 speci c cold-regulated DEGs (associated with T10), 41 DEGs were enriched in the KEGG pathways, and most of them were classi ed to following pathways: 'Metabolic pathway', 'Carbon metabolism', 'Citrate cycle (TCA)', 'RNA transport', and 'Lysosome'. Interestingly in T20 and T40, '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 25 , and the a chrysomelid beetle, Galeruca daurica 16 . Transcriptome analysis revealed that the expression of 'Glycolysis' and 'TCA cycle pathways' are up-regulated in a similar manner as the braconid wasp, Aphidius colemani, when exposed to low temperatures 29 . 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 22 . 'Fatty acid synthase' and 'Fatty acid metabolism', are two of the main pathways of RIFA when exposed to high temperature. Expression of fatty acids as hydrophobic agents allows insects to avoid water loss in warmer regions of the globe 30 . At high temperatures, Gomphocerus sibiricus are known to increase their levels of oleic acid, linoleic acid, linolenic acid and glycerin, and phenomenon can suppress mortality due to excessive evaporation of body moisture 31 . 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 28 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 32 . The synthesis and metabolism of amino acids are necessary to produce a signi cant number of amino acids, which make available the raw materials necessary for the synthesis of heat-resistant proteins 28 .
In this study, two cuticular protein unigenes were identi ed from 203 co-regulated DEGs under cold temperature stress (Table S9). Cuticular protein gene expression has been observed in studying other insects studies such as beetles, moths, planthoppers, and stick insects when exposed to cold temperature stress 16,22,25,33 . Although the physiological role of cuticular proteins in insect cold hardiness has not yet been identi ed, it seems insect cuticle may play an important role in insects when coping with low temperature 16,19,22,34,35 .
According GO analysis (Fig. 6A), 'Antioxidant activity' was enriched at low temperatures. Suggesting that is might contribute to RIFA ability to resist oxidative stress damage at low temperature 16 , or their potential for cell preservation via antioxidant defense when in challenged by environmental complexity 36 . In addition, one cytochrome P450 was identi ed that up-regulated exclusively under low temperature (Table  S9). Meanwhile NADH dehydrogenase subunit 1 was up-regulated at all treated temperature (Fig. 5C, Table S8). These two proteins are main enzymes at antioxidant activity pathway 25 . In 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 22 .
Heat shock proteins (HSPs) are another important protein that insects use as critical physiological products when under abiotic stress conditions 28 . From earlier studies it was believed that Hsp is associated with biological cold and heat resistance 37,38 . HSPs are molecular chaperones, which 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 39,40 . In this study, we identi ed one speci c Hsp70 that up-regulated about a 12-fold change when RIFA was exposure 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 con rmed that the expression of Hsp genes can be up-regulated by cold and heat stimulus 39,41 . To assist the resistance to temperature stress, the Hsp60 gene expression in Stegobium paniceum signi cantly increases under highand low-temperature stress 42 . Three Hsp90 and four Hsp70 were up-regulated by cold stress and were differentially expressed at desert beetle, Microdera punctipennis 21 . The differences in Hsp, insect species, sex of organism, and intensity of temperature are important factors related to Hsp expression level in insects 16,43 .
In conclusion, we compared the transcriptomes of S. invicta under high-and low-temperature stresses using RNA-Seq technology based on the high-throughput sequencing. Comparative transcriptome analysis identi ed many genes, and a large number of changes were discovered in the metabolic pathways through the GO and KEGG enrichment analysis. Our data will facilitate further molecular investigations and genomic research. Many novel relationships between high-and low-temperature and signi cantly up-regulated genes were identi ed in this study (Table S7- , and an arti cial diet described by Dussutour and Simpson (2008) 44 . 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, 30, and 40℃ for 24 h. Untreated control ants were incubated at 30℃. 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℃ for subsequent experiments.

RNA extraction and RT-qPCR
RNAs samples were extracted from the whole body of S. invicta adults using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. After RNA extraction, it was resuspended in nuclease-free water and quanti ed using a spectrophotometer (NanoDrop, Thermo Scienti c, Wilmington, DE, USA). cDNA was then synthesized from RNA (1 µg) using RT PreMix (Intron Biotechnology, Seoul, Korea) containing oligo dT primer according to the manufacturer's instruction. All quantitative PCRs (qPCRs) in this study were determined using a real time PCR machine (CFX Connect Real-Time PCR Detection System, Bio-Rad, Hercules, CA, USA) and iQ SYBR Green Supermix (Bio-Rad, Hercules, CA, USA) according to the guidelines provided by the manufacturer. The reaction mixture (20 µL) contained 10 µL of iQ SYBR Green Supermix, 1 µL of cDNA template (100 ng), 1 µL each of forward and reverse primers (Table S13), and 7 µL nuclease free water. RT-qPCR cycling began with 95°C heat treatment for 10 min followed by 40 cycles of denaturation at 94°C for 30 s, annealing at 52°C for 30 s, and extension at 72°C for 20 s. The expression level of Ef1_β as a reference gene was used to normalize target gene expression levels 45 under different treatments. PCR products were assessed by melting curve analysis. Quantitative analysis was performed using comparative CT (2 −ΔΔCT ) method 46 .

Illumina sequencing
To obtain short-read RNA sequences, Illumina sequencing was performed at Macrogen (Seoul, Korea). Each library was constructed from 1 µg total RNA from the whole body of 5 individuals (not pooled) of S. invicta adults per treatment using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, USA) and sequenced using HiSeq 4000 System (Illumina, San Diego, USA) with 101 bp pair end read (Table S1).

De Novo Assembly
Illumina short reads were quality-ltered and adapter-trimmed using Trimmomatic v0.38 (http://www.usadellab.org/cms/?page=trimmomatic). FastQC v0.11.7 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to check data quality before and after trimming. After the removal of low-quality reads, an Illumina-based de novo transcriptome assembly was performed using Trinity version trinity rnaseq r20140717, bowtie 1.1.2 47 . Trimmed reads for every sample were merged into one le 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 48 . 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 ltered and clustered into the nonredundant transcripts using CD-HIT version 4.6 (http://weizhongli-lab.org/cd-hit) 49 . These transcripts were de ned as 'unigenes' which are used for predicting the ORFs (Open Reading Frames), annotating against several known sequence databases, and analyzing differentially expressed genes (DEGs). ORF prediction for unigenes was performed using TransDecoder version 3.0.1 (https://github.com/TransDecoder/TransDecoder/wiki) 50 to identify candidate coding regions within transcript sequences. After extracting ORFs that were at least 100 amino acids long, TransDecoder predicted the likely coding regions. Trimmed reads for each sample were aligned to the assembled reference using Bowtie program. For the differentially expressed gene analysis, the abundances of unigenes across samples were estimated into read count as an expression measure by RSEM algorithm

Differential Gene Expression Analysis
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. In order to reduce systematic bias, we estimated the size factors from the count data and applied Relative Log Expression (RLE) normalization with 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. Signi cant unigene results were analyzed as Up and Down-regulated count by log 2 FC ≥ 5, ≤ -5 and ≥ 10, distribution of expression level between two groups was plotted as Volcano plot and simple bar plots.

Quantitative RT-PCR validation
The twenty genes in response to cold treatment (T10) were chosen for validation using qRT-PCR. To do that, S. invicta adults were incubated at 10° and 30°C for 24 hr in two separate groups that included 10 ants. RNA extraction and cDNA synthesis were performed according to 'RNA extraction and RT-qPCR' section. Speci c primers were designed using Primer Quest tool (www.idtdna.com) (Table S13). Expression level of Ef1_β was used as a reference gene and to normalize target gene expression levels under different treatments 45 . PCR products were assessed by melting curve analysis. Quantitative analysis was performed using comparative CT (2 −ΔΔCT ) method 46 . Finally, the data were compared according ratio of FPKM and ratio of mRNA expression levels for all selected genes. Declarations