Different bacterial and viral pathogens trigger distinct immune responses in a globally invasive ant

Invasive species populations periodically collapse from high to low abundance, sometimes even to extinction. Pathogens and the burden they place on invader immune systems have been hypothesised as a mechanism for these collapses. We examined the association of the bacterial pathogen (Pseudomonas spp.) and the viral community with immune gene expression in the globally invasive Argentine ant (Linepithema humile (Mayr)). RNA-seq analysis found evidence for 17 different viruses in Argentine ants from New Zealand, including three bacteriophages with one (Pseudomonas phage PS-1) likely to be attacking the bacterial host. Pathogen loads and prevalence varied immensely. Transcriptomic data showed that immune gene expression was consistent with respect to the viral classification of negative-sense, positive-sense and double-stranded RNA viruses. Genes that were the most strongly associated with the positive-sense RNA viruses such as the Linepithema humile virus 1 (LHUV-1) and the Deformed wing virus (DWV) were peptide recognition proteins assigned to the Toll and Imd pathways. We then used principal components analysis and regression modelling to determine how RT-qPCR derived immune gene expression levels were associated with viral and bacterial loads. Argentine ants mounted a substantial immune response to both Pseudomonas and LHUV-1 infections, involving almost all immune pathways. Other viruses including DWV and the Kashmir bee virus appeared to have much less immunological influence. Different pathogens were associated with varying immunological responses, which we hypothesize to interact with and influence the invasion dynamics of this species.

(DWV), and Kashmir bee virus (KBV). These "honey bee" viruses have also been observed to infect and replicate in other insects 15,16 , though Argentine ants appear to have higher viral loads when interacting with honey bees 14 . At any one time colonies of ants or bees are frequently co-infected with a range of viral and other pathogens referred to as the 'pathosphere' [17][18][19][20] . Pathogens may interact to influence other pathogens within their hosts. In addition to the viruses, Argentine ants have a diverse bacterial community that spatially varies both internationally and even within countries like New Zealand 21 . Very little is known about the influence of these viruses, bacterial species, or the microbial pathosphere on Argentine ant physiology.
The immune response of any ant species to viral infection is poorly understood. Much more is known about how honey bees respond to viruses and other pathogens. McMenamin et al. 22 reviewed honey bee and insect immune pathways. They described how the genome of honey bees encodes the major components of insect immune pathways, with a focus on genes influenced by viral infections. Gene expression patterns altered by viruses include those in the immune pathways: Jak/STAT (Janus kinase/Signal Transducer and Activator of Transcription), RNAi (RNA interference), Toll, Imd, and JNK (c-Jun N-terminal kinase). For example, viruses are known to increase expression of genes such as hopscotch (Jak/STAT pathway), dicer and argonaute-2 (RNAi) in honey bees 23 . Infection of honey bees by different viruses from different viral families results in a pattern of immune gene expression with some common genes, but also in many distinct expression patterns 23 . The immune response to viral challenges in ants and bees also appears to be influenced by the presence of other pathogens, including bacteria 12,18,24 . Recent work in Argentine ants has indicated there may be a core set of immune genes that are influenced by infection with pathogens such as bacteria. These core immune genes include many in the Toll and Imd pathways 25 . Different pathogens could compete for host resources or may even be synergistic in their influence on the ant or bee physiology.
Our goal in this study was to examine the association of variation in viral community or pathosphere with immune gene expression in Argentine ants. We first used an RNA-seq approach to examine expression patterns of immune genes in Argentine ant communities with varying viral species and infection loads. We then focused on how four viruses are associated with the expression of four different immune pathways, represented by nine immune genes. Because the immune response to viruses may be influenced by other microbial pathogens, we also included in our analysis a putative bacterial pathogen (Pseudomonas spp.) that we previously identified as common in Argentine ants 21 . Pseudomonas bacteria are major pathogens in a wide range of both plants and animals 26,27 . Overall, we wanted to understand how different viruses are associated with the immune response in these globally invasive ants. Do all the viral pathogens have a similar influence on the ant's immune physiology? How are different immune pathways within the ants influenced by different viruses and bacteria?

Results
Virus and immune gene analysis via RNA-seq. We sampled ant populations for the RNA-seq analysis only from sites in Northland (Fig. 1a), so as to limit the influence of extraneous variables on the ant microbiota and gene expression. Seventeen viruses were putatively found in the RNA-seq analysis, including evidence for three different bacteriophage viruses. Six viruses were found in all samples analysed (i.  toti-like virus 1, for example, showed a ≤ 10-fold level of variability in infection between samples. In contrast, the two viruses showing the most variability in viral abundance between RNA-seq samples were Linepithema humile C-virus 1 with a 12,032-fold change between highest and lowest infection, and the Deformed Wing Virus with a 5,968-fold change between highest and lowest infection. Two other viruses that we assayed in the RT-qPCR analysis were also amongst the top five most variable virus loads between samples (Linepithema humile virus 1 with 464-fold variation, and Kashmir Bee Virus with 62-fold variation). Our previous work has suggested that virus loads may change in ants when they are in the presence of honey bees 14 . The RNA-seq sampling did not provide sufficient data to examine this question statistically as, for example, only two of eight samples tested positive for KBV in sites without honey bees. The available data indicated that the range of infection loads was similar for many viruses in sites with and without honey bees. Some viruses, however, such as DWV and LHUV-1 had viral loads that were orders of magnitude higher when in the presence of honey bees ( Table 1).
The PLS-regression analysis aimed at finding associations between the detected viruses and known Argentine ant immune genes. Hierarchical clustering of viruses indicated that immune gene expression was generally consistent with respect to the viral classification of negative-sense, positive-sense and double-stranded RNA viruses (Fig. 2). The patterns of gene expression from infection were most distinct in association with infection by the two negative-sense RNA viruses LhuRV1 (Linepithema humile rhabdo-like virus 1) and LhuBLV1 (Linepithema humile bunya-like virus 1), while most positive-sense as well as dsRNA viruses appeared to elicit similar immune responses. LhuCV1 (Linepithema humile C-virus 1), a positive-sense RNA virus, did not show a comparable pattern to any other virus in terms of its association with immune gene expression. A large number of genes across all immune pathways were positively correlated with negative-sense RNA viruses. Interestingly, the same set of genes appeared down-regulated in response to dsRNA and most positive-sense viruses (see the top half of the heatmap). Surprisingly, genes that function within the RNAi and JaK/STAT pathways, which have been previously associated with anti-viral immune defence in bees 24 were often negatively associated with positive-sense and dsRNA viruses, and appeared to be mostly activated by negative-sense viruses. In contrast, genes that were the most strongly associated with positive-sense RNA viruses were peptide recognition proteins assigned to the Toll and Imd pathways (e.g. beta-1,3-glucan-binding and peptidoglycan recognition proteins) and displayed negative associations with negative single-stranded RNA viruses.

Virus, bacteria and immune gene analysis by RT-qPCR.
Our RT-qPCR analysis examined viruses, bacteria and immune gene expression in samples of 30 ants sampled from around New Zealand. Of the viruses, the most prevalent to least prevalent were LHUV-1 (49 of 52 samples), followed by DWV (43/52), BQCV (32/52), and KBV (30/52). Comparison between species of viral infection loads is difficult due to PCR amplification constraints, though of all viruses LHUV-1 and DWV had an order of magnitude higher average level of infection than KBV and BQCV. BQCV appeared to have the lowest rate of infection or load within ants. Within viruses we observed a substantial amount of variation between samples. For DWV, there was a 276.3 million-fold increase in viral loads from lowest to most infected sample. Levels of variation were similarly high for LHUV-1 (a 16.7 million-fold change) and KBV (an 8.4 million-fold change). The bacteria Pseudomonas spp. was positive in all www.nature.com/scientificreports www.nature.com/scientificreports/ samples. It ranged in relative expression from 0.131-103.390 (a 788-fold level of variation). The comparison of sites with and without honeybees indicated significantly higher viral loads for BQCV, DWV and KBV in apiaries ( Fig. 1b-d). In contrast, LHUV-1 viral loads in ants were orders of magnitude higher outside of apiaries, and there was no difference in Pseudomonas spp. loads (Fig. 1e,f). There were several significant correlations between viral infection loads. For example, LHUV-1 was negatively correlated with KBV (r s = −0.47, P < 0.001, Fig. 3). BQCV was positively correlated with KBV (r s = 0.47, P < 0.001, Fig. 3).
We observed much less variation in the expression of the immune genes. The nine genes showed a 3-to 12-fold level of variation between the minimum and maximum values. Spearman rank correlations indicated that the expression of any immune gene was typically highly correlated with other immune genes (Fig. 3). This result was expected for genes within pathways. For example, dicer and Argonaut are both within the RNAi pathway and their expression in ants was highly correlated (r s = 0.84, P < 0.001, Fig. 3). Correlations between different immune genes were significantly positive and negative, likely indicating pathways have a mixture of activation of and inhibition of expression.
Our principal components analysis (PCA) showed that the five original, Box-Cox transformed (but unstandardized) pathogens were almost linearly independent: each linearly independent component had a single pathogen coefficient that explained the vast majority (≥93.1%; or ≥0.965 2 ) of the component's variation (Table 2). The Heatmap showing the regression matrix between viral loads for 11 viruses and gene expression for 90 immune genes obtained after Partial Least Squares regression analysis. Virus types clustered together (top dendrogram, colour-coded in the top bar), suggesting that specific sets of immune genes are specifically associated with positive and negative single-stranded, and double stranded RNA viruses. We included immune genes associated with the Toll, Imd, RNAi, Jak/STAT and JNK immune pathways, colour-coded in the left-hand bar). The colour key from blue to yellow to red indicates the strength of the regression between each virus and immune gene.
www.nature.com/scientificreports www.nature.com/scientificreports/ ordering of the principal components from most to least variance gave an effective ordering of the variance in pathogen loads: the bacteria Pseudomonas spp., LHUV-1, DWV, KBV and BQCV. We then estimated a regression model for each immune gene on the five pathogen principal components (Table 3). These regression models suggest that the Pseudomonas spp. and LHUV-1 were important for the majority of the immune genes. Both DWV and KBV each significantly associated with only one of the nine immune genes. BQCV was not observed to be significantly associated with the expression of any immune gene. Note that the coefficients need ceteris paribus interpretations: the results for each immune gene are in response to a unit change in each pathogen principal component while holding all other variables constant, and should not be interpreted in isolation. The principal components are linearly independent by construction though, so there is no correlation between any of the predictors but the signs of the coefficients in each component are arbitrary (as in any PCA). Interestingly, however, www.nature.com/scientificreports www.nature.com/scientificreports/ the signs of the coefficients for the majority of the immune genes in Pseudomonas spp. were the opposite to those of LHUV-1 and that contrast is meaningful.

Single ant analysis of viruses and immune genes by RT-qPCR.
In order to further understand viral loads and the Argentine ant immune system, we also examined virus and immune gene expression in individual Argentine ants. All ants were collected from the same site, colony and time. The viruses showed an extremely high level of variation, with a 160-fold level of variation for LHUV-1 and 30-fold for DWV (Table 4). In contrast the immune genes defensin-2 (DEF2), tyrosine-protein kinase hopscotch (HOP), and argonaute-2 (AGO) showed relatively lower levels of variation in expression between ants.

Discussion
We found that invasive Argentine ants in New Zealand have a diverse and highly variable pathosphere. Viral infections varied between sites by several orders of magnitude and were associated with a variable immune response. Of the viruses examined in detail via RT-qPCR, the Deformed wing virus (DWV) and the Linepithema humile virus 1 (LHUV-1) were highly variable with some ants showing a 276 million-fold difference in load between sites. The immune response appeared to be greatest to the bacteria Pseudomonas spp. and LHUV-1 infections, while much reduced for BQCV.
The field of eco-immunity is relatively new and still developing. Cornet et al. 28 highlight how host immunity is at the core of major theories for invasion biology. Changes in pathogen communities can result in a reallocation of resources from immune defences to beneficial traits associated with increased reproductive output or invasive potential. The potential for such a reallocation makes sense if different pathogens and pathogen communities "Pseudomonas" "LHUV-1" "DWV" "KBV" "BQCV"   Table 3. Linear regression results showing the estimated coefficients (standard errors) of the relationships between the immune genes and the five principal components, shown here as the five pathogens from the analysis presented in Table 2 www.nature.com/scientificreports www.nature.com/scientificreports/ invoke varying costs for their hosts 29 . The presence of honey bees appeared to be associated with a substantial alteration of the viral community within ants (although we acknowledge it is difficult to be conclusive in this result as ant-bee interactions were only observed in northern sites). We observed variation in infection and immune responses throughout the sampling sites. In other ants, different viral pathogens have previously been demonstrated to incur different phenotypic effects. For example, SINV-1 significantly reduced fire ant (Solenopsis invicta) queen weight but not reproductive output, while SINV-2 reduced reproductive output but not weight 10 . Each Solenopsis invicta virus appears to induce a different immune response in their ant host 10 . Transcriptomic data analysis indicates that different viruses similarly appeared to induce a different response in Argentine ants, with responses largely grouped by viral classification. Interestingly, PCA analysis indicates that different viruses of similar types appeared to be associated with different immune responses. This observed variation in immune response to different viruses might be induced by a range of different mechanisms, which include some virus families having molecular machinery that inhibits certain antiviral immune responses 30 .
Microbial species including gut microbiota have been hypothesized to influence immune responses, pathogen communities, and ultimately their host's invasion success 28 . Perhaps the most extreme example of bacteria influencing viral infections in insects is the intracellular bacteria Wolbachia inhibiting infection of mosquitoes from the virus causing dengue fever 31 . An extreme example of viruses influencing bacterial infections are phages, which have even been described as mutualists for the hosts of pathogenic bacteria 32 . Our study provides the first evidence of phages in Argentine ants. Pseudomonas is a Gram-negative bacterial pathogen with a wide host range that includes plants, insects and animals 26,27 . It can be a virulent bacterial pathogen of insects with an LD 50 of less than 10 bacterial cells 26 . Pseudomonas has been found in Argentine ants 21,33,34 , although it is yet to be confirmed as an ant pathogen. Gram-negative bacteria have been associated with an increased immune response along the Imd pathway 24 , which we observed, though this bacterium also appeared to influence most of the immune genes that we quantified. Interestingly, the sign of the coefficients for the response of the immune genes to Pseudomonas was typically the opposite for that to LHUV-1 (Table 3). Pathogen loads within ant samples were positively correlated for some species, and negatively correlated for others. Different microbial pathogens are clearly associated with different immune responses, which may be related to variation in replication strategies and virulence among different classes of virus, or among different viruses. For instance, in mosquitoes, the Toll pathway is efficient against Flavivirus, but not Alphavirus, both positive single-stranded RNA viruses 35 . Given that different virus species within the same class can illicit such different immune responses, it is easy to conceive that different sets of genes will be affected by different viral classes. Viljakainen et al. 25 similarly found viruses and bacteria to jointly influence the immune response of Argentine ants, via both up-and down-regulation of different genes. Overall, it appears the model of a 'pathosphere' 18-20 seems appropriate here for Argentine ants. The pathogen community likely interacts and jointly influences the host's immune response, with the phages perhaps influencing some bacteria. While our models explained a high proportion of the variation in immune gene expression, other microbial taxa undoubtedly influenced the health and immune response of these ants. Argentine ants have a diverse bacterial community 21 and other viruses 12 that we did not quantify, or that are still unknown to science.
Currently little is known about the influence of these viruses on ants. For the viruses and bacteria that we analysed, our results suggest that these pathogens are associated with a strong immune response, especially LHUV-1 and Pseudomonas spp. A comparatively much-reduced response was observed in association with DWV, KBV and BQCV infection levels. Samples of ants showed high levels of variation in the viral loads. For example, in the composite samples of 30 ants in the RT-qPCR analysis, there was a 16.7 million-fold increase in LHUV-1 loads from lowest (non-zero) to most infected sample. The overall levels of variation are consistent with variation in viral loads observed in honey bees 36,37 , bumble bees 36 , and wasps 38 . Our analysis of within-site variation for  Table 4. An analysis of patterns of viral load and immune gene expression in single ants, all of which were sampled from the same colony at the same time. Virus loads were extremely variable between ants, relative to the variation observed in immune gene expression. The viruses were Linepithema humile virus 1 (LHUV-1) and the Deformed wing virus (DWV). The immune genes examined were defensin-2 (DEF2), tyrosine-protein kinase hopscotch (HOP), and argonaute-2 (AGO).
www.nature.com/scientificreports www.nature.com/scientificreports/ single ants showed a 28-fold variation for DWV, and a 161-fold difference between the highest and lowest level of LHUV-1 infection. These levels of variation between ants are of some concern for sampling and attempts to understand the influence of viruses on ant populations. If high viral titres are only observed in a small percentage of ants within a population then a single ant with, for example, a ~200-fold higher viral load in some individuals will skew estimates of viral loads.
Variation in observed immune gene expression was much less than for viruses. The samples of 30 pooled ants showed immune genes varied by 3-to 12-fold. For single ants, only a 3-to 5-fold level of variation was observed for the three genes we examined (defensin-2, tyrosine-protein kinase hopscotch, and argonaute-2). It is noteworthy that the patterns of expression for nearly all the immune genes were correlated both within and between immune pathways. We expected that dicer and argonaute expression would be highly correlated because they are within the RNAi pathway. Surprisingly, nearly all the other genes associated with the other immune pathways examined were also significantly and highly correlated with dicer and argonaute (Fig. 3). Some genes were positively and some negatively correlated, as has previously been observed in pathogen infection in Argentine and other ants 10,25 . The inter-relatedness of the immune system response, bacterial and viral infections supports an eco-immunological framework 28,39 and the hypothesis that the microbial community within Argentine ants can exert an immunological cost. The extent of that immunological cost in relation to variable pathogen pressure, and whether or not it can induce population collapses of this invader 40 , or promote variation in its ecological success, remains to be determined.

Methods
Ant collection and RNA extraction. Ants were collected using an aspirator from eighteen locations across New Zealand (Fig. 1a). In locations to the north of New Zealand, Argentine ants were collected from honey bee (Apis mellifera) apiary sites where they were observed infesting hives and attacking bees. All ants were collected alive, immediately snap frozen in liquid nitrogen, and then stored at −70 °C until RNA extraction. RNA for the RNA-seq analysis was extracted using Direct-Zol RNA Microprep kit for the RNA-seq experiment, and Direct-Zol RNA Miniprep for the RT-qPCR experiment (Zymo Research, CA, USA). Briefly, ants were homogenized in 300, 600 and 100 µL of TRIzol Reagent (Life Technologies/ThermoFisher Scientific) for the Microprep, Miniprep for samples of pooled ants and Miniprep for single ant extractions, respectively. Pooled samples consisted of 30 ants from single sites. Homogenates were clarified by centrifugation and RNA purified following manufacturer's instructions. For the RNA-seq samples, final RNA concentrations and integrity were determined using a BioAnalyzer (Agilent technologies, USA) and the RNA stored at −70 °C until samples were selected and shipped for sequencing. For the RT-qPCR samples, we determined RNA concentrations using a NanoDrop One (ThermoFisher Scientific) and likewise stored the samples at −70 °C until processed.
Virus and immune gene analysis via RNA-seq. Twelve samples were selected for RNA sequencing, based on a preliminary RT-qPCR analysis of virus infection using RNA extracted from 30 ants per sample. We included samples collected in sites with and without honeybee hives, and selected samples that spanned a wide range of viral loads. RNA samples were stored and shipped in RNAStable (Biomatrica, USA) for sequencing by Annoroad (China). Libraries were prepared using the NEBNext Ultra RNA library Prep Kit for Illumina according to the manufacturer's instructions and 150 bp paired-end reads were sequenced on a single lane using Illumina HiSeq platform. Data pre-processing included the removal of read pairs when either one read contained adapter contamination of more than 5 nucleotides, contained more than 10% of ambiguous nucleotide calls or contained more than 50% of low-quality base calls (Phred score <20, i.e. 99% accuracy).
Clean paired-end reads were aligned to the Argentine ant reference genome 41 with HISAT 2.1.0 with default parameters 42 to produce sample-specific BAM files. We then fed the BAM output into StringTie 1.3.4 43 to generate GTF files, using the -e argument to restrict the assembly to transcripts matching the Argentine ant RefSeq genome records (GCF_000217595.1). We generated a raw transcript counts matrix at the gene level using the authors' prepDE.py script. We imported the raw counts matrix into R as a DGEList object using the edgeR 3.22.3 package 44 , together with one of the GTF files created by StringTie as a source of information on each gene. In order to filter the raw counts matrix for low expressed transcripts we first computed counts per million (CPM), then discarded transcripts that were not expressed in at least 3 samples at more than 1 CPM. We then applied a TMM normalisation on the filtered raw count data. Expression data in CPM for 90 annotated genes likely associated with immune pathways were extracted and transformed into FPKM using transcript lengths, and inputted as a matrix for PLS-regression analysis. We selected the immune genes based on the KEGG Pathway database as well as searching the Argentine ant RefSeq genome records for genes known to be associated with immunity in closely related organisms. Of all 90 genes included in the analysis, 40 (44%) were in the Toll immune pathway, 30 (33%) were in the Imd pathway, 11 (12%) in the RNAi pathway, 6 (7%) in the JNK pathway, and 3 (3%) in the Jak/STAT pathway.
Reads that did not align to the Argentine ant genome in the HISAT step were de novo assembled using Trinity 2.3.2 with default parameters. We quantified assembled transcript expression within Trinity using eXpress 1.5.1 (Roberts and Pachter 2012), yielding a TMM-normalised TPM matrix at the gene level. We then aligned the assembled transcripts to various reference databases in order to screen for and quantify viral transcripts: we used BLAST 2.2.3 45 to run BLASTx searches of the NCBI viral protein database as well as of a more recently published Argentine ant novel viruses database 12 , using a e-value cut-off of 10 −5 . We discarded hits that were less than 100 bp long and returned less than 90% identity with the target sequence. From the filtered BLASTx output, we selected a single best hit per assembled transcript based on the highest bitscore. In order to compute viral loads, we first summed TPM values for all genes belonging to each identified virus. Then, we used the Argentine ant library size to normalise viral loads to host tissue among samples. Log-transformed TMM-normalised transcripts per million (TPM), normalised to Argentine ant read counts, were stored as a matrix for later analysis.
www.nature.com/scientificreports www.nature.com/scientificreports/ In order to explore associations between Argentine ant immune response and viral pressure, we used the immune gene expression and viral load matrices as input for a Partial Least Square analysis (PLS) in regression mode with eight components, using mixOmics 46 . We plotted a heatmap of similarity matrix showing the association between viral loads and immune gene expression, based on the two components that were significant. All remaining statistical analyses were performed in R 3.5.1 47 . Virus, bacteria and immune gene analysis by RT-qPCR. cDNA was generated from 500 ng RNA/ sample using qScript XLT cDNA SuperMix (QuantaBio, MA, USA) in 20 µL reactions. The cDNA was diluted to a final volume of 212 µL with water and 8 µL/sample loaded into each well of the MicroAmp Fast Optical 96-Well Reaction Plate (Applied Biosystems/ThermoFisher Scientific) for RT-qPCR. Each reaction also contained: 10 µL PowerUp SYBR Green Master Mix (Applied Biosystems/ThermoFisher Scientific), and forward and reverse primers at final concentrations of 300 nM. Plates were analysed on a QuantStudio 7 Flex Real-Time PCR System (Applied Biosystems) using fast PCR cycling conditions (50 °C, 2 min; 95 °C, 2 min; 40 cycles of 95 °C, 1 s, 60 °C, 30 s). Cycle threshold (Ct) values were used to calculate target (Ct target ) levels relative to the reference (Ct reference ) genes using the equation (2^(−Ct target ))/(2^(−Ct reference )). We used the R package lmPerm (permutation ANOVA for linear models; Wheeler & Torchiano 48 ) to examine for differences in relative virus loads between populations of Argentine ants collected near bee hives and from areas without hives. The targets of the primers used for the above analysis, presence, identity, and, for selected species, viral replication, were confirmed via Sanger sequencing.
Performance Analytics 49 was used within the R statistical environment to plot the data, following Box-Cox transformations on each variable in the car package 50 to help reduce the influence of extreme values. Spearman correlation coefficients were calculated between immune gene expression and viral/bacterial loads, with a sequential Bonferroni adjustment for the p-values. We used principal components analysis to represent the unscaled variation in the five Box-Cox transformed pathogen loads as five linearly independent (latent) components, then linear regression models were used to relate the response of each immune gene to the five independent pathogen components.