Transcriptomic basis and evolution of ant nurse-larval social regulatory interactions

Development is often strongly regulated by interactions among close relatives, but the underlying molecular mechanisms are largely unknown. In eusocial insects, interactions between caregiving worker nurses and larvae regulate larval development and resultant adult phenotypes. Here, we begin to characterize the social interactome regulating ant larval development by collecting and sequencing the transcriptomes of interacting nurses and larvae across time. We find that the majority of nurse and larval transcriptomes exhibit parallel expression dynamics across larval development. We leverage this widespread nurse-larva gene co-expression to infer putative social gene regulatory networks acting between nurses and larvae. Genes with the strongest inferred social effects tend to be peripheral elements of within-tissue regulatory networks and are often known to encode secreted proteins. This includes interesting candidates such as the nurse-expressed giant-lens, which may influence larval epidermal growth factor signaling, a pathway known to influence various aspects of insect development. Finally, we find that genes with the strongest signatures of social regulation tend to experience relaxed selective constraint and are evolutionarily young. Overall, our study provides a first glimpse into the molecular and evolutionary features of the social mechanisms that regulate all aspects of social life. Author Summary Social interactions are fundamental to all forms of life, from single-celled bacteria to complex plants and animals. Despite their obvious importance, little is known about the molecular causes and consequences of social interactions. In this paper, we study the molecular basis of nurse-larva social interactions that regulate larval development in the pharaoh ant Monomorium pharaonis. We infer the effects of social interactions on gene expression from samples of nurses and larvae collected in the act of interaction across a developmental time series. Gene expression appears to be closely tied to these interactions, such that we can identify genes expressed in nurses with putative regulatory effects on larval gene expression. Genes which we infer to have strong social regulatory effects tend to have weak regulatory effects within individuals, and highly social genes tend to experience relatively weaker natural selection in comparison to less social genes. This study represents a novel approach and foundation upon which future studies at the intersection of genetics, behavior, and evolution can build.

Given that we observed transcriptome-wide patterns consistent with nurse-larva transcriptional co-154 regulation across larval development, we next identified the genes that might be driving these patterns 155 (see S3 Fig). We performed differential expression analysis to identify genes that varied in larval 156 expression according to larval developmental stage, as well as genes that varied in nurse expression 157 according to the developmental stage of larvae they fed. We identified 8125 differentially expressed 158 genes (DEGs) in larvae (78% of 10446 total genes). We identified 2057 and 1408 DEGs in stage-specific 159 nurse heads and abdomens, respectively, compared to 599 and 520 DEGs in random nurse heads and 160 abdomens, respectively. We removed genes differentially expressed in both stage-specific and random 161 nurses (N = 272 DEGs in heads, N = 140 DEGs in abdomens), which might differ among our colony 162 replicates due to random colony-specific effects that were not consistently associated with social 163 regulation of larval development. After this removal, we retained the top 1000 DEGs, sorted by P-value, 164 for each sample type other than random nurses (larvae, stage-specific nurse heads, stage-specific nurse 165 abdomens) for social gene regulatory network reconstruction, reasoning that these genes were the most 166 likely to be involved in the regulation of larval development. 167 168

Reconstruction of social gene regulatory networks 169
To infer putative gene-by-gene social regulatory relationships between nurses and larvae, we 170 reconstructed gene regulatory networks acting within and between nurses and larvae (S3 Fig). The output 171 of regulatory network reconstruction is a matrix of connection strengths, which indicate the regulatory 172 effect (positive or negative) one gene has on another, separated according to the tissue the gene is 173 expressed in. To identify the most highly connected (i.e. centrally located, upstream) genes of regulatory 174 networks, we calculated within-tissue connectivity and social connectivity by averaging the strength of 175 connections across each connection a gene made, differentiating between within-tissue (nurse-nurse or 176 larva-larva) and social connections (nurse-larva) (Fig 1B). On average, within-tissue connectivity was 177 higher than social connectivity (Wilcoxon test; P < 0.001 in all tissues), and within-tissue connectivity 178 was negatively correlated with social connectivity in each tissue (S4 Fig). The top enriched gene ontology terms based on social connectivity in nurses were entirely dominated by metabolism (S3,S4 Tables; see  180   also S5 Table for the top 20 genes by nurse social connectivity). 181 182 Secreted proteins and social gene regulation 183 While based on our data it is not possible to distinguish between genes that code for protein products that 184 are actually exchanged between nurses and larvae versus genes that affect behavior or physiology within 185 organisms (Fig 1A), proteins that are known to be cellularly secreted represent promising candidates for 186 the social regulation of larval development [40]. We downloaded the list of proteins that are known to be 187 cellularly secreted from FlyBase [52] and used a previously-generated orthology map to identify ant 188 orthologs of secreted proteins [40]. Genes coding for proteins with orthologs that are cellularly secreted in 189 Drosophila melanogaster had higher social connectivity than genes coding for non-secreted orthologs in 190 nurse heads (Fig 3A; Wilcoxon test; P = 0.025), though not for nurse abdomens (P = 0.067). 191 For the most part, we have focused on broad patterns of nurse-larva gene coregulation. In this 192 paragraph, we will highlight the potential social role of one of the genes with the highest social 193 connectivity within nurse heads, giant-lens (S6 Table; giant-lens is the 7 th highest gene coding for 194 secreted proteins by social connectivity in nurse heads). Giant-lens is an inhibitor of epidermal growth 195 factor receptor (EGFR) signaling [53], and giant-lens expression in nurse heads was negatively associated 196 with the expression of the homolog of eps8, human EGFR substrate 8 in larvae, most prominently seen in 197 the spike in nurse giant-lens expression accompanied by a drop in larval eps8 expression at the end of 198 larval development (Fig 3B). Giant-lens was also used in regulatory network reconstruction in larvae (i.e. It is important to note that these patterns were not seen for 203 all genes in the EGFR pathway, and the results presented here cannot be taken as concrete evidence of 204 EGFR regulation via social processes. Nonetheless, the mechanism illustrated here represents a tangible 205 example of how nurse-larva interactions could function at the molecular level. 206 207 Molecular evolution of social gene regulatory networks 208 To investigate the selective pressures shaping social regulatory networks, we used population genomic 209 data from 22 resequenced M. pharaonis workers, using one sequenced M. chinense worker as an outgroup 210 [48]. Using polymorphism and divergence data, we estimated gene-specific values of selective constraint, 211 which represents the intensity of purifying selection that genes experience [54]. To identify genes 212 disproportionately recruited to the core of social regulatory networks, we calculated "sociality index" as 213 the difference between social connectivity and within-tissue connectivity for each gene. Sociality index 214 was negatively correlated to selective constraint due to a positive correlation between within-tissue 215 connectivity and constraint and a negative correlation between social connectivity and constraint ( Fig 4A-216 C). Additionally, genes differed in sociality index according to their estimated evolutionary age, with 217 ancient genes exhibiting lower sociality indices than genes in younger age categories ( Fig 4D). Finally, 218 while evolutionary age and evolutionary rate appear to be somewhat empirically confounded [55], 219 selective constraint and evolutionary age were each independently associated with sociality index, based 220 on a model including both variables as well as tissue (GLM; LRT; evolutionary age: c 2 = 21.536, P <

224
In organisms with extended offspring care, developmental programs are controlled in part by socially-225 acting gene regulatory networks that operate between caregivers and developing offspring [14,42]. In this 226 study, we sequenced the transcriptomes of ant nurses and larvae as they interacted across larval 227 development to assess the effects of social interactions on gene expression dynamics. We found that large 228 sets of genes (i.e. modules) expressed in ant larvae and their caregiving adult nurses show correlated changes in expression across development (Fig 2). The majority of nurse and larval transcriptomes was 230 represented in these correlated modules, suggesting that the tight phenotypic co-regulation characterizing 231 nurse-larva interactions over the course of larval development is also reflected at the molecular level. 232 To characterize the overall network and evolutionary patterns of genes involved in nurse-larva 233 interactions, we reverse engineered nurse-larva gene regulatory networks and calculated the "social 234 connectivity" for each gene, defined as the sum of inferred social regulatory effects on all genes 235 expressed in social partners. We found that genes with high social connectivity tended to have low 236 within-individual connectivity (S4 Fig; where within-individual connectivity is defined as the sum of 237 inferred regulatory effects acting within a given tissue). Nurse-expressed genes with higher sociality 238 indices (i.e disproportionately higher social connectivity than within-individual connectivity) tended to be 239 evolutionarily young and rapidly evolving due to relaxed selective constraint (Fig 4). Genes with high 240 social connectivity were enriched for a number of Gene Ontology (GO) categories associated with 241 metabolism (S3,S4 Tables), consistent with the idea that molecular pathways associated with metabolism 242 are involved in the expression of social behavior [56,57]. Previously, many of the proteins found to be 243 widely present in social insect trophallactic fluid transferred from nurses to larvae were involved in sugar 244 metabolism (e.g. Glucose Dehydrogenase, several types of sugar processing proteins) [15]. Along the 245 same lines, many of the genes with with high social connectivity in our study are also annotated with 246 terms associated with sugar metabolism (S5 Table; e.g. Glycerol-3-phosphate dehydrogenase, Glucose 247 dehydrogenase FAD quinone, Pyruvate dehydrogenase). Finally, we found that genes encoding for 248 orthologs of cellularly-secreted proteins in Drosophila melanogaster (possibly important for intercellular 249 signaling) tended to exhibit higher levels of social connectivity than their non-secreted counterparts (Fig  250   3A). 251 One gene that stands out in terms of being cellularly secreted and exhibiting a relatively high 252 social connectivity is giant-lens, which inhibits EGFR signaling [ [61]. Further experimental work is necessary to ascertain 256 whether giant-lens is actually orally secreted by nurses and transferred to larvae, but gene expression 257 dynamics are consistent with the social transfer of giant-lens from nurses to larvae, followed by the 258 inhibition of EGFR signaling at the end of larval development in worker-destined larvae ( Fig 3B). In general, the suite of genic characteristics commonly associated with worker-biased genes 271 (rapidly evolving, evolutionarily young, loosely connected) are all consistent with relaxed selection acting 272 on genes associated with workers [49]. Here, we show that within the worker caste, genes that appear to 273 be functionally involved in the expression of social behavior (i.e. nursing) experience relaxed selective 274 constraint relative to genes important for within-worker processes. Therefore, the combination of kin 275 selection as well rapid evolution thought to be characteristic of social traits [25] likely act in concert to 276 shape the labile evolutionary patterns commonly associated with worker-biased genes. Finally, it has also 277 been suggested that plastic phenotypes such as caste recruit genes which were evolving under relaxed 278 selection prior to the evolution of such plastic phenotypes [70-72]. Our results could also be consistent 279 with this hypothesis, though the population genomic patterns we observe show that relaxed selective 280 constraint is ongoing. 281 In this study, we sought to reconstruct regulatory networks acting between nurses and larvae, 282 beginning with the assumption that nurse gene expression changes as a function of the larval stage fed. 283 This is more likely to be the case when nurses are specialized on feeding particular larval stages. 284 According to a previous study, about 50% of feeding events are performed by specialists (though note 285 specialization is likely a continuous trait, and the 50% figure is the result of a binomial test) [40]. 286 Therefore, we expect our stage-specific nurse samples to comprise about 50% specialists. We also expect 287 random nurse samples to contain 50% specialist nurses, but, crucially, the specialists should be relatively 288 evenly divided among larval stages since random nurses were collected regardless of which larval stage 289 they were observed feeding. Because our stage-specific nurse samples did not consist of 100% specialists, 290 we expect that the signal of nurse-larva co-expression in our analysis is effectively diluted. In order to 291 maximize the signal of nurse-larval co-expression dynamics, future studies would ideally focus entirely 292 on specialists, as well as on tissues such as brains and the specific exocrine glands [73] known to be 293 important for social behavior and communication. Despite these limitations, we were still able to observe 294 transcriptomic signatures consistent with the social regulation of larval development. 295 296

297
In this study, we uncovered putative transcriptomic signatures of social regulation and identified distinct 298 evolutionary features of genes that underlie "social physiology", the communication between individuals 299 that regulates division of labor within social insect colonies [74,75]. Because we simultaneously collected 300 nurses and larvae over a time series of interactions, we were able to elucidate the putative molecular 301 underpinnings of nurse-larval social interactions. This is a promising approach that could be readily 302 extended to study the molecular underpinnings of all forms of social regulation in social insect colonies, 303 including regulation of foraging, regulation of reproduction, etc.. Furthermore, by adapting the 304 methodology presented here (i.e. simultaneous collection over the course of interactions followed by 305 sequencing), the molecular mechanisms and evolutionary features of genes underlying a diverse array of 306 social interactions, including courtship behavior, dominance hierarchy formation, and regulation of biofilm production could all be investigated. Overall, this study provides a foundation upon which future 308 research can build to elucidate the genetic underpinnings and evolution of interacting phenotypes. 309 310

311
This study builds on previous work investigating genomic signatures of kin selection in which we 312 characterized transcriptomic profiles from adult queens and workers, as well as queen-and worker-313 destined larvae [48]. While stage-specific nurses were used in the previous analysis, the knowledge of the 314 developmental stage of larvae they fed was not, as they were simply treated as adult workers. This study 315 also complements the past dataset with new data from random nurses, which were collected concurrently 316 with previous samples. 317

Study Design 319
To construct experimental colonies, we began by creating a homogenous mixture of approximately fifteen 320 large source colonies of the ant Monomorium pharaonis. From this mixture, we created thirty total 321 replicate experimental colonies of approximately equal sizes (~300-400 workers, ~300-400 larvae). We 322 removed queens from ½ the study colonies to promote the production of reproductive-destined larvae. We pre-assigned colonies to one of five larval developmental stages (labeled L1-L5, where L1 333 and L2 refer to 1st-instar and 2nd-instar larvae and L3, L4, and L5 refer to small, medium, and large 3rd-334 instar larvae [77]). We identified larval stage through a combination of hair morphology and body size. 335 L1 larvae are nearly hairless, L2 larvae have straight hairs and are twice the length of L1 larvae, and L3-336 L5 larvae have dense, branched hairs [78]. We separated 3rd-instar larvae into three separate stages based 337 on body size [77] because the vast majority of larval growth occurs during these stages. We sampled 338 individuals (larvae as well as nurses) across larval development time: beginning at the L1 stage, we 339 sampled colonies assigned to each subsequent stage at intervals of 3-4 days, by the time the youngest 340 larvae in colonies lacking queens were of the assigned developmental stage (note that in colonies lacking 341 queens, no new eggs are laid so the age class of the youngest individuals progressively ages). We sampled 342 each colony once, according to the developmental stage we had previously assigned the colony (e.g. for 343 colonies that we labeled 'L4', we waited until it was time to sample L4 larvae and nurses and sampled 344 individuals from that colony at that time). From each colony, we sampled stage-specific nurses and 345 worker-destined larvae, as well as random nurses from colonies with queens and reproductive-destined 346 larvae from colonies without queens (starting at the L2 stage, because at L1 caste cannot be distinguished 347 For each time point in each assigned colony, we collected stage-specific nurses, nurses feeding 351 larvae of the specified developmental stage (L1, L2, etc). Concurrently, we collected random nurses, 352 nurses we observed feeding a larva of any developmental stage. Rather than paint-marking nurses, we 353 collected them with forceps as soon as we saw them feeding larvae. We collected random nurses as soon 354 as we observed them feeding a larva of any developmental stage in the course of visually scanning the 355 colony. We did not make an attempt to systematically collect nurses from different areas of the nest but 356 did so haphazardly, such that the distribution of larval stages fed resembled overall colony demography. 357 Nurses feed L1 and L2 larvae exclusively via trophallaxis (i.e. liquid exchange of fluid), while nurses feed L3-L5 larvae both via trophallaxis and by placing solid food in larval mouthparts [79]. To get a 359 representative sample of all types of nurses, we did not distinguish between nurses feeding liquid and 360 solid food, though all L3-L5 samples contained a mixture of the two. After collecting nurses, we 361 anaesthetized the colony using carbon dioxide and collected larvae of the specified developmental stage.  We used the package EdgeR [82] to construct models including larval developmental stage and replicate 385 and performed differential expression analysis for each sample type separately. We retained genes 386 differentially expressed according to a nominal P-value of less than 0.05 (i.e. no false discovery 387 correction), as the purpose of this step was simply to identify genes that could be involved in interactions 388 that shape larval development (rather than spurious interactions arising from replicate-specific effects). 389 See S1 Dataset for a list of all stage-specific nurse and larval differentially expressed genes. As a side note, a version of GENIE3 that was developed for time series data, dynGENIE3 [85], 407 does exist. However, we opted to utilize the original GENIE3 algorithm because we reasoned that the 408 temporal spacing of developmental stages was likely too sparse for regulatory network reconstruction to 409 incorporate time (note also that the co-expression algorithm we used, STEM, was explicitly designed for 410 short time series such as ours). While our method therefore does not explicitly incorporate temporal 411 dynamics, we purposefully biased our results to emphasize larval development over differences between 412 replicates by only utilizing genes differentially expressed across larval development (or based on larval 413 stage fed in the case of nurses). 414 We repeated the entire regulatory reconstruction reconstruction process 1000 times and averaged 415 pairwise connection strengths across runs, as the algorithm is non-deterministic. To capture the total 416 effect of each gene on the transcriptome dynamics within tissues, we averaged the regulatory effects each 417 gene had on all other 999 genes expressed in the same tissue ("within-individual connectivity"). 418 Similarly, to capture the effect each gene had on the transcriptome of social partners, we averaged 419 regulatory effects each gene had on the 1000 genes expressed in social partners ( "social connectivity").