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).