Experimental increase in fecundity causes upregulation of fecundity and body maintenance genes in the fat body of ant queens

Temnothorax rugatulus is a small ant with colonies of a few hundred workers and one to several queens. Two queen morphs can occur and we only used colonies of the common larger queen morph (19). These ants reside in rock crevices in forests throughout Western North America. We collected 105 T. rugatulus colonies in the Chiricahua Mountains, Arizona in August 2015 (Table S1). In our laboratory, colonies were transferred to artificial nest boxes and kept at 22°C and 12h light / 12h dark in a climate chamber. For the dietary-restriction experiment, we limited the queens’ access to workers, as these might buffer food restrictions imposed on the queen. The queen was isolated with five workers to ensure food provisioning in the upper part of an artificial nestsite (queenright part, QR), while the reminder of the colony inhabited the lower section (queenless, QL). Both parts were separated by a metal grid, allowing the exchange of volatiles, but not food (Fig. S1). The QR-parts of the dietary-restriction treatment were provided with two cricket legs and a droplet of honey every 2nd week (N = 32 colonies, Fig. S1), while the QR-parts of the control (N = 30 colonies) and the QL-parts of both treatments the received the same amount of food twice weekly. At each feeding session, all nests were opened and any remaining food was removed. Food was replaced with fresh food at every session for the QL-parts and the QR-part control, but only at every 4th time for the QR-part of dietary-restriction treatment. Thus, these ants had food available only a quarter of the time, whereas all others had continuous access to food. All colonies had continuous access to water. Worker survival was monitored weekly. Queens reduced egg-laying over time in both treatments, as eggs laid at the beginning developed into larvae. In order to increase the likelihood to detect an treatment effect, all eggs and young larvae were removed from the QR parts at week eight (Figure S2). The experiment ended after 13 weeks. All eggs were counted. All queens were dissected. For the egg-removal experiment, 44 polygynous colonies were used to create 58 experimental colonies. 14 colonies were split and colony fragments were allocated to different treatments. We standardized the number of queen, workers, and larvae to 2, 50 and 12 respectively, and removed all eggs. In the egg-removal treatment (N = 29 colonies) all eggs were removed once per week, while in the control treatment (N = 39 colonies) eggs were just moved around with forceps. Colonies were anesthetized with CO2 to remove or simulate egg removal. Queen survival was recorded weekly. The experiment was performed over six weeks. All queens were dissected two weeks after the experiments’ end. Ovaries and fat bodies of eight queens per treatment were dissected on ice (N = 32). The fat body was individually homogenized in 50μL TRIZOL (Invitrogen) and stored at -20°C. RNA was extracted using the RNeasy mini kit (Qiagen) with a preceding chloroform step. Library preparation and sequencing of 100bp paired-end reads on an Illumina HiSeq 2000/2500 was conducted at BGI Hong Kong. The ovaries of all remaining queens were dissected and photographed for fertility measurements (Leica DFC425 20x; LAS version 4.5). Ovary length in the dietary-restriction experiment was analysed by using a Wilcoxon test. We used generalized-linear models with a poisson distribution (link function = log) to investigate the effect of treatment on the number of white eggs in the ovaries and in the colony as dependent variables. For the egg-removal experiment, fecundity differences were analysed with a linear-mixed model with ovary length (in mm) as dependant variable, and a generalized linear-mixed model with a poisson distribution (link function = log) with the number of white eggs in the ovaries as dependent variable. Experimental fragment ID and colony ID were added as random factors. For both experiments, we separately analysed queen survival by running survival models, with treatment as explanatory variable. As all queens were independent in the dietary-restriction experiment, we used the R package survival (), while we used the package coxme() for the analysis of the egg removal data by adding colony ID as random factor. The statistical analyses were conducted in R v. 3.0.2 (R Development Core Team 2008). For the transcriptome analyses, raw reads from all 32 samples were trimmed with Trimmomatic-v0.36 (20), quality checked using FastQC-v0.11.5 (21). Paired reads were de-novo assembled using Trinity v.2.4.8 (22), resulting in 328,731 transcripts. For annotation, we conducted a BlastX homology search (23) against the non-redundant invertebrate protein database (June 2018) with an E-value cut-off of E-05. Read count estimates per transcript and sample were obtained using RSEM-v1.3.0 (24) with Bowtie2 aligner for each experiment separately. To eliminate low read counts likely representing noise, we removed transcripts with less than 10 reads in less than four samples (25). The differential expression analyses were performed with R package Deseq2-v1.2.10 (26) (contrast function) by comparing treatment to control for each experiment. We added colonyID as random factor in the egg-removal treatment as some samples were dependent. Nucleotide sequences were translated into amino-acid sequences with Transdecoder-v5.5.0 (22), before conducting a gene ontology (GO) annotation using InterProScan-v5.34-73.0 (27). We performed GO term enrichment analyses based on subsets of DEGs using the R package TopGo-v-3.6 (28), with the “weight01” algorithm. For each DEG, we extracted the geneID from the BlastX results to retrieve additional GO and biological functions from the Uniprot database (www.uniprot.org) with Homo sapiens, Mus musculus, and Drosophila melanogaster as query organisms using an in-house python script (Supplement: maintenance_test.R). Thereafter, we searched for terms associated with fecundity (fecund, fertile, meiosis, meiotic, zygote, reproductive, reproduction, embryo, pregnancy, mating, foetal, sexual, brood, egg, ovule, ovary, ovarian), body maintenance (Toll, response to oxidative stress, apoptosis, TOR, tumour repressor, transposable element, response to UV damages, DNA repair, stress response, aging, autophagy, cellular homeostasis), epigenetics (chromatin, histone), fatty acid metabolism (fatty), and immunity (immune). χ²-tests were used to contrast the frequency of DEGs with these functions between treatment and controls. We conducted these additional analyses to obtain insights into putative functions of DEGs in fecundity, longevity (body maintenance, immunity), food processing (fatty acid metabolism) and gene regulation (epigenetics).

红足盘腹蚁（Temnothorax rugatulus）是一种小型蚂蚁，其蚁群包含数百只工蚁与1至数只蚁后。该物种存在两种蚁后形态，本研究仅采用常见的大体型蚁后形态的蚁群（19）。红足盘腹蚁栖息于北美西部各地森林的岩石缝隙中。本研究于2015年8月在亚利桑那州奇里卡瓦山脉采集了105个该蚁种的蚁群（补充表S1）。
在实验室环境中，我们将蚁群转移至人工巢箱，并置于温度22℃、光周期12小时光照/12小时黑暗的气候培养箱中饲养。为开展饮食限制实验，我们限制了蚁后与工蚁的接触，因为工蚁可能会缓冲施加给蚁后的食物限制压力。我们将1只蚁后与5只工蚁一同隔离在人工巢的上部区域（有蚁后区域，QR），以保障该区域的食物供给，而蚁群的其余个体则栖息于下部区域（无蚁后区域，QL）。两个区域通过金属网格分隔，仅允许挥发性物质交换，无法传递食物（补充图S1）。
饮食限制处理组的有蚁后区域每两周投喂2根蟋蟀腿和1滴蜂蜜（共32个蚁群，补充图S1）；而对照组的有蚁后区域（共30个蚁群）以及两个处理组的无蚁后区域则每周投喂两次等量食物。每次投喂前，我们会打开所有蚁巢并清除残留食物。无蚁后区域与对照组有蚁后区域每次投喂均更换新鲜食物，而饮食限制处理组的有蚁后区域仅每四次投喂才更换一次新鲜食物。因此，该组蚁群仅在四分之一的时间内可获得食物，其余组群则可持续获取食物。所有蚁群均可持续获得饮水。
研究人员每周监测工蚁的存活情况。两组处理中的蚁后产卵量均随时间推移而下降，因为早期产下的卵已发育为幼虫。为提升检测处理效应的概率，我们在第8周时移除了有蚁后区域内的所有卵与幼龄幼虫（补充图S2）。该实验于13周后结束，对所有蚁卵进行计数，并对所有蚁后进行解剖。
在蚁卵移除实验中，我们使用44个多蚁后蚁群构建了58个实验蚁群：其中14个蚁群被拆分，蚁群片段被分配至不同处理组。我们将每个实验蚁群的蚁后、工蚁与幼虫数量分别标准化为2、50和12，并移除所有蚁卵。蚁卵移除处理组（共29个蚁群）每周移除所有蚁卵，而对照组（共39个蚁群）仅用镊子移动蚁卵以模拟移除操作。使用二氧化碳对蚁群进行麻醉，以便移除或模拟移除蚁卵。每周记录蚁后的存活情况，该实验周期为6周。实验结束两周后，对所有蚁后进行解剖。
每个处理组取8只蚁后的卵巢与脂肪体，在冰上进行解剖（共32个样本）。将脂肪体分别置于50μL TRIzol（Invitrogen公司）中匀浆，并保存于-20℃环境。使用RNeasy小型试剂盒（Qiagen公司）提取RNA，步骤前需加入氯仿处理。文库构建与100bp双端测序在Illumina HiSeq 2000/2500平台上完成，测序工作由香港华大基因（BGI Hong Kong）执行。对其余所有蚁后的卵巢进行解剖并拍照，以测量其繁殖能力（使用Leica DFC425 20倍显微镜；LAS V4.5软件）。
饮食限制实验中的卵巢长度数据采用Wilcoxon检验进行分析。我们采用服从泊松分布的广义线性模型（连接函数为对数函数），以蚁巢内与卵巢内的白色卵粒数量为因变量，分析处理方式对其的影响。在蚁卵移除实验中，我们以卵巢长度（单位：mm）为因变量，采用线性混合模型分析繁殖能力差异；以卵巢内白色卵粒数量为因变量，采用服从泊松分布的广义线性混合模型（连接函数为对数函数）进行分析。将实验片段ID与蚁群ID作为随机效应加入模型。
两项实验均单独通过生存模型分析蚁后存活情况，以处理方式作为解释变量。由于饮食限制实验中的蚁后均为独立个体，我们使用R语言survival包进行分析；而蚁卵移除实验中，我们通过加入蚁群ID作为随机效应，使用coxme()包完成数据分析。所有统计分析均在R v3.0.2环境中完成（R开发核心团队，2008年）。
转录组分析方面，我们使用Trimmomatic-v0.36（20）对32个样本的原始测序reads进行修剪，使用FastQC-v0.11.5（21）进行质量检测。利用Trinity v2.4.8（22）对双端reads进行从头组装，最终获得328731条转录本。注释阶段，我们以E值≤1e-5为阈值，针对非冗余无脊椎动物蛋白质数据库（2018年6月版）进行BlastX同源性搜索（23）。针对每个实验，我们分别使用RSEM-v1.3.0（24）结合Bowtie2比对工具，估算每个样本中每条转录本的读段计数。为去除可能代表测序噪声的低读段计数转录本，我们移除了在少于4个样本中读段数小于10的转录本（25）。
差异表达分析使用R语言DESeq2-v1.2.10包（26）的对比函数完成，针对每个实验比较处理组与对照组。由于蚁卵移除实验中部分样本存在相关性，我们将colonyID作为随机效应加入模型。使用TransDecoder-v5.5.0（22）将核苷酸序列翻译为氨基酸序列，随后使用InterProScan-v5.34-73.0（27）进行基因本体（Gene Ontology, GO）注释。我们使用R语言TopGo-v3.6包（28）结合“weight01”算法，针对差异表达基因（Differentially Expressed Genes, DEGs）子集进行GO功能富集分析。
针对每个差异表达基因，我们从BlastX搜索结果中提取geneID，并通过自研Python脚本（补充材料：maintenance_test.R）在UniProt数据库（www.uniprot.org）中检索其额外GO注释与生物学功能，检索的参考物种包括智人（Homo sapiens）、小家鼠（Mus musculus）与黑腹果蝇（Drosophila melanogaster）。随后，我们检索与以下功能相关的术语：繁殖能力（fecund、fertile、meiosis、meiotic、zygote、reproductive、reproduction、embryo、pregnancy、mating、foetal、sexual、brood、egg、ovule、ovary、ovarian）、机体维持（Toll、氧化应激反应、细胞凋亡、TOR、肿瘤抑制因子、转座因子、UV损伤响应、DNA修复、应激反应、衰老、自噬、细胞稳态）、表观遗传学（chromatin、histone）、脂肪酸代谢（fatty）以及免疫（immune）。我们采用卡方检验（χ²检验）对比处理组与对照组中具有上述功能的差异表达基因的频率。
我们开展上述额外分析，以探究差异表达基因在繁殖能力、寿命（机体维持与免疫）、食物代谢（脂肪酸代谢）以及基因调控（表观遗传学）中的潜在功能。