Results
Effects of social chromosome genotype on gene expression. To analyze the effect of the supergene-bearing social chromosome’s genotype on gene expression patterns in S. invictagynes, we first utilized a generalized linear model accommodating variance due to the experimental factors (colony identity, natal social form, tissue type, genotype, and covariation between tissue type and genotype) in our RNA-seq data (multifactorial-model). We found that the different social chromosome genotypes are associated with many highly significant differences in gene expression independent of tissue type, and that presence of the inversion-based Sb haplotype leads to a relative increase in gene expression at most of the significantly affected loci (FDR < 0.01; Figures 2A,B and S1). This pattern is particularly prominent in the SB/SB versus SB/Sbcomparison, where 96% of the differentially expressed genes (DEGs) are upregulated in heterozygous gynes relative to homozygous SB/SBgynes (Figures 2A and S1). It is less pronounced in the SB/SBversus Sb/Sb comparison (72% of DEGs are upregulated inSb/Sb individuals), although there are more than double the number of DEGs in this comparison (as expected given the greater difference in Sb copy number in the latter versus the former comparison; Figures 2B and S1). These data are consistent with a scenario in which the Sb haplotype drives widespread upregulation of gene expression, but such a mechanism would appear to be moderated in the absence of a copy of the SB haplotype.
We next tested for organ- and genotype-specific effects of the supergene using pairwise comparisons exclusively between polygyne-reared individuals. We observed identical numbers of significant DEGs in the brains and ovaries in comparisons of SB/SB versus SB/Sbgynes (42 DEGs in each; FDR < 0.01; Figures 2C and S1). More than double this number was observed in both organs in the SB/SBversus Sb/Sb comparisons, again as expected given the greater difference in Sb copy number in the latter comparison, with the highest number of such DEGs occurring in the ovaries (Figures 2C and S1). Minimal functional enrichment for GO terms was observed in any of the comparisons, but this could be due in part to the lack of functional information for many genes in the S. invicta genome (Supplementary Data 3).
We observed a strong bias in the localization of DEGs to the supergene region in the comparisons of different social chromosome genotypes (Figure 2D). More than a quarter of these DEGs occur within the supergene region despite this region accounting for only 4% of mapped gene space in our assembly (X2-test; P < 0.01). Moreover, we observed distinct patterns of allelic dominance with respect to social chromosome genotype for DEGs within the supergene region as compared to the rest of the genome (Figure S3, Watson’s Two-Sample Test of Homogeneity, P < 0.05). These results suggest that cis -regulatory divergence is an important outcome of sequence divergence between the Sb and SB supergene haplotypes. However, we also observed widespread trans -acting effects of Sb presence, with between 49% and 78% of DEGs attributable to social chromosome genotype located outside of the supergene region (Figures 2D and S2). This is noteworthy because, while the Sb and SB supergene regions rarely undergo recombination and consequently bear many fixed genetic differences, the remainder of the genome recombines freely and shows minimal sequence differentiation between males bearing or lacking Sb (Yan et al., 2020).
Effects of natal social form on gene expression. Monogyne and polygyne colonies of S. invicta differ in many ways other than colony queen number, for example adult worker to brood ratios, dispersion and connectivity of nests, and population densities (Tschinkel 2006). The cuticular semiochemical profiles and reproductive capacity of queens differ between the two forms, presumably driven largely by their social chromosome genotypes (C. J. DeHeer, 2002; Eliyahu, Ross, Haight, Keller, & Liebig, 2011), yet indirect genetic effects stemming from the different supergene compositions of the colony worker force in the two forms also potentially affect the phenotypes of developing gyne offspring, independently of gyne genotype. For example, substantial effects of the social developmental environment on gene expression were previously observed in workers of S. invictacolonies (Wang, Ross, & Keller, 2008). Thus, we expected to see differences in gene expression between gynes reared in these very different types of societies even if they possess the same social chromosome genotype. Surprisingly, we found that no genes were significantly differentially expressed solely in response to social form of origin when assessed across tissues by a multifactorial-model (Figure 3A). Similarly, pairwise comparisons of specific tissues fromSB/SB gynes from different natal colony types (monogyne or polygyne) revealed only 22 social form-dependent DEGs (FDR < 0.01), and these were confined to brains (Figure 3B). However, we did observe an apparent synergistic effect between natal colony social form and social chromosome genotype on differential expression in the ovaries. We observed 43 DEGs between ovaries of SB/SB andSB/Sb individuals when both were reared in polygyne colonies but this increased to 205 DEGs when polygyne-reared SB/Sb individuals were compared to monogyne-reared SB/SB individuals (Figure 3C).
Allele-specific expression in the supergene region. The high proportion of genes upregulated in SB/Sb heterozygotes relative to SB/SB homozygotes that occur within the supergene region can be explained by two potential mechanisms operating withinSB/Sb individuals: 1) elevated expression of the Sballeles relative to alleles from the SB homologous region, or 2) a relatively balanced increase in the expression of both the Sband SB alleles. Elevated expression of Sb alleles would be consistent with differences in SB and Sb chromatin structure, mutations in supergene regulatory elements that disproportionately influence Sb -linked alleles in cis , or the presence of Sb -specific paralogs (Fontana et al., 2019). A balanced increase in SB and Sb alleles would be consistent with Sb -induced activation of regulatory elements that can interact with both alleles of a gene, or with Sb -induced chromatin remodeling that extends to both alleles. In order to distinguish between these mechanisms, we evaluated allele-specific expression (ASE) in heterozygotes at genes within the supergene region by means of diagnostic SNPs (SB versus Sb , derived from our RNA-seq data). Our analysis identified 28 genes (of the 222 with sufficient diagnostic SNP coverage) exhibiting significant ASE (FDR < 0.01, Figure 4A, Supplementary Data 2). Notably, we observed similar frequencies of significant ASE for the SB - andSb -linked alleles, as reported previously in S. invictausing a more limited experimental design (Wang et al., 2013). However, of the nine genes that exhibit significantly elevated expression inSB/Sb relative to SB/SB gynes, six also exhibit significant Sb ASE (Figure 4B), providing support for the importance of elevated expression of the Sb alleles relative to alleles from the SB homologous region. The other three of these genes exhibit differential expression without significant ASE (Figure 4B), serving as examples of relatively balanced upregulation ofSB - and Sb -linked alleles in heterozygotes.
Our analysis of ASE within the supergene region of the genome allowed us to investigate putative SB dosage compensation in heterozygous individuals. Genes subject to such dosage compensation are expected to lack significant differential expression between SB/SB andSB/Sb gynes (FDR > 0.01) but to exhibit significantSB allele-specific overexpression (P < 0.01). We identified twelve such genes (Supplementary Data 2), raising the prospect that SB dosage compensation operates at these loci, possibly due to diminished functionality of the homologous Sbvariants.
Heterogeneity of differential expression in the supergene region. In order to gain insight into the regulatory topology of the supergene region, we assessed whether DEGs and genes with ASE are distributed uniformly throughout the supergene region or located in discernible clusters. Plots of gene expression pattern by social chromosome position qualitatively illustrate that genes exhibiting either type of differential expression are widespread yet spatially clustered according to effect within the supergene region (Figure 5, S5, Supplementary Data 2). To statistically assess whether DEGs are randomly distributed throughout the region, we employed Wald-Wolfowitz runs tests (Wald & Wolfowitz, 1940) on binary representations of significant DEGs, specifying whether gene expression is elevated in the presence or absence of the Sbhaplotype and ordering DEGs based on relative linear position in the supergene region (Supplementary Data 2). For DEGs between SB/SBand SB/Sb gynes, only one SB/SB- upregulated gene maps to the supergene and, unsurprisingly, the test is nonsignificant (P > 0.05). In contrast, among DEGs between SB/SB andSb/Sb individuals, we found that genes with significantly elevated expression in the presence or absence of the Sbhaplotype, respectively, are not distributed randomly throughout the supergene region (Wald-Wolfowitz runs test, P = 0.019). This heterogeneous pattern persists when the analyses are repeated using a recently published, improved genome assembly (Figure S5, Supplementary Data 4) (Yan et al., 2020).
In order to illustrate the contrasting patterns of gene regulation observed within the fire ant supergene, we assigned genes to two exemplary regions. The first contains loci that exhibit a high degree ofSB allele-specific expression coupled with SB/SB -elevated expression in the SB/SB versus Sb/Sb comparisons (region 1); the second is characterized by Sb allele-specific expression coupled with elevated expression in heterozygotes and Sbhomozygotes relative to SB homozygotes (region 2). The appearance of these distinct regions shows that supergene haplotype divergence has generated complex, non-uniform changes in the transcriptional landscape.
Candidate genes associated with polygyny. To glean functional insights into the molecular regulation of colony social form in S. invicta , we grouped strongly supported DEGs by patterns of differential expression into four categories: 1) those with expression differences between polygyne-reared SB/SB and SB/Sb gynes that were brain-specific , 2) those with differences between gynes of these same classes that were ovary-specific , 3) those with differences between SB/SB and Sb/Sb gynes but notSB/SB and SB/Sb gynes, independent of natal social form or organ type (Sb/Sb-specific) , and 4) those with differences in all comparisons of Sb- supergene presence/absence, independent of natal social form or organ type (constitutive ). We hypothesized that brain-specific DEGs would exist because of the dramatically different dispersal, mating, and colony-founding behaviors of gynes bearing or lacking the supergene. Similarly, we hypothesized that ovary-specific DEGs would exist because of the dramatic differences in rates of ovary development and associated physiological processes (e.g., fat and storage protein deposition) of young queens of different supergene status. Given the effective lethality of the Sb/Sbgenotype in gynes (Hallar et al., 2007), we expected that genes differentially expressed solely in the SB/SB-Sb/Sb comparisons might give some indications as to the molecular causes of the recessive lethal effect. Finally, the genes always differentially expressed when comparing individuals with and without the Sb haplotype, given their ubiquitous differential expression, are prime candidates for drivers of the multifaceted Sb- mediated polygyne phenotype. A list of genes in each category is presented in Table 1 (a complete description of expression statistics for each gene and a table of genes that appeared as DEGs in multiple analyses but do not fit into one of our four candidate gene categories is provided in Supplementary Data 2).