Mating and blood-feeding induce transcriptome changes in the spermathecae of the yellow fever mosquito Aedes aegypti

Aedes aegypti mosquitoes are the primary vectors of numerous viruses that impact human health. As manipulation of reproduction has been proposed to suppress mosquito populations, elucidation of biological processes that enable males and females to successfully reproduce is necessary. One essential process is female sperm storage in specialized structures called spermathecae. Aedes aegypti females typically mate once, requiring them to maintain sperm viably to fertilize eggs they lay over their lifetime. Spermathecal gene products are required for Drosophila sperm storage and sperm viability, and a spermathecal-derived heme peroxidase is required for long-term Anopheles gambiae fertility. Products of the Ae. aegypti spermathecae, and their response to mating, are largely unknown. Further, although female blood-feeding is essential for anautogenous mosquito reproduction, the transcriptional response to blood-ingestion remains undefined in any reproductive tissue. We conducted an RNAseq analysis of spermathecae from unfed virgins, mated only, and mated and blood-fed females at 6, 24, and 72 h post-mating and identified significant differentially expressed genes in each group at each timepoint. A blood-meal following mating induced a greater transcriptional response in the spermathecae than mating alone. This study provides the first view of elicited mRNA changes in the spermathecae by a blood-meal in mated females.

Mating alters the transcript profile in the spermathecae. To investigate the global transcriptomic response to mating and blood-intake in spermathecal tissue, we first analyzed the combined post-mating timepoints. This approach reduced complexity, simplified the dimensionality of the data, and reduced the false discovery rate, enabling us to identify candidate transcripts in NBF and BF females. We found 275 transcripts that were significantly differentially expressed in spermathecae from both NBF and BF females compared to virgins ( Fig. 2; Figure S2A-B; File S3)-107 up-and 168 down-regulated ( Fig. 2A). We observed a higher number of differentially expressed transcripts in spermathecae from BF females (152 transcripts: 54 up-and 98 downregulated; Fig. 2A, right panel), compared to NBF females (121 transcripts: 51 up-and 70 down-regulated; Fig. 2A, left panel). A similar number of transcripts were uniquely up-regulated in spermathecae from NBF and BF females; 10 and 13 transcripts, respectively (Fig. 2B). However, we observed a significant difference in the number of down-regulated genes in spermathecae from BF females compared to NBF females (DF = 1; χ 2 = 33.5; P < 0.0001; Fig. 2B). These results indicate that mating and blood-feeding each induce a transcriptional response in the spermathecae. In hematophagous insects such as Ae. aegypti, blood intake causes significant physiological and metabolic constraints 50,51 . Thus, while spermathecae from BF females and NBF females overlap in their physiological and metabolic processes for sperm maintenance during storage, stress caused by a blood-meal might independently effect a transcriptional response in the spermathecae.
Regardless of female feeding status, sperm require energetic nourishment, protection from oxidative stress and a suitable pH environment within storage 20 . Our results show a total of 40 transcripts that are up-regulated in response to mating in the spermathecae of NBF and BF females ( Fig. 2B; File S3). The most significant up-regulated genes in spermathecae from NBF and BF females correspond to proteins with roles in ion transport, such as bumetanide-sensitive sodium coupled cation-chloride cotransporter (AAEL009888) and sodium-dependent phosphate transporter protein (AAEL001656) ( Fig. 2A; File S3). Ion transport genes are also regulated by mating in the spermathecae of social insects 29 and in the D. melanogaster seminal receptacle 32 , a second sperm storage organ in this species, indicating the importance of ion balance in insect sperm storage. As insect cells have mechanisms to control pH through ion exchange (e.g. Na, H, K, etc.) 52,53 , the up-regulation of these genes might be related to the ionic balance and/or pH of the storage organs. Moreover, we found other transcripts potentially related to the regulation of the spermathecae's osmotic environment, such as putative calcium ion binding transcripts troponin C (AAEL012114), EF-hand domain-containing protein (AAEL006006), and several genes whose products relate to transmembrane transport (AAEL027912, AAEL010883, AAEL024256, AAEL000488, AAEL012041, AAEL007979, AAEL012041, AAEL014311, AAEL020162) (File S3).
Sperm maintenance within storage requires energy provisioning substrates such as polysaccharides and lipoproteins 20 . Although energy substrates for sperm maintenance have not been identified in mosquitoes 54 , it is reasonable to hypothesize that energy for sperm nourishment may be provided by the spermathecal cells 21 . We found two transcripts in the spermathecae from both NBF and BF females that encode proteins that potentially function in sperm nourishment and energy metabolism: the oxysterol-binding protein related to lipid transport (AAEL018354) and phospholipase A1 (AAEL006982). In spermathecae from BF females, we found genes related  www.nature.com/scientificreports/ to energy provisioning such as trehalase (AAEL009658) and transcripts related to amino-acid metabolism such as asparagine synthetase (AAEL015631) and cysteine dioxygenase (AAEL026357). Moreover, during sperm storage, the imbalance of reactive oxygen species and antioxidants can cause oxidative stress 55 . Antioxidant genes and proteins have been reported in A. mellifera 39 and An. gambiae 23 spermathecae, and the D. melanogaster 32 seminal receptacle. We found three transcripts with potential roles in oxidative stress in the BF and NBF female groups: egg-antigen (AAEL001094), the synaptic vesicle membrane protein VAT-1 homolog-like (AAEL019604), and a predicted oxidative stress-induced growth inhibitor (AAEL002835) ( Fig. 2A; File S3). The egg-antigen belongs to the heat-shock protein family. Heat shock proteins have roles in reducing oxidative stress in other insects, including Drososphila 56 and the honey bee 57 .
Differential expression analysis also identified several down-regulated transcripts in NBF and BF females. A total of 55 genes are significantly down-regulated in spermathecae from both NBF and BF females compared to virgins (Fig. 2B). The most significant down-regulated transcripts correspond to cytochrome P450 49a1 (AAEL008638), a zinc finger protein (AAEL004057), sodium/hydrogen exchanger 9B2 (AAEL011109) and nardilysin (AAEL010073) (File S3). Post-mating down-regulation of P450 genes-which encode proteins with several functions including oxidation catalysis of exogenous substances-occurs in the D. melanogaster sperm storage organs 32 . Regulation of P450 transcripts might minimize oxidative stress and suppress an immune response to protect sperm cells within storage 32,41,58 . Temporal transcriptional profile dynamics of the spermathecae in response to mating and blood-feeding. Physiological and structural changes in female reproductive tract tissue and/or stored sperm occurs in insects during and soon after mating ends [59][60][61][62][63][64][65] . In D. melanogaster, oviductal tissue remodeling 61 , uterine conformation 65,66 , and increased innervation of oviduct musculature 63 are detectable in mated females during, or by 6 h after mating. In An. gambiae, permanent structural changes of the female reproductive tract are detectable at 8 h post-mating 62,67 . Aedes aegypti sperm undergo a coat modification within the spermathecae detectable at 4 h and complete by 24 h post-mating 64 . Mating responsive genes in the female reproductive tract potentially mediate post-mating tissue or sperm modification and/or coordinate activities with seminal proteins that are required for some of these changes 59,60,63 . Transcriptome studies of Drosophila female reproductive tract tissues show that peak differential gene expression occurs at 6 h post-mating 68,69 . In An. gambiae, the number of genes that experience a twofold or greater change in expression increases at each consecutive post-mating time point assessed (2, 6, and 24 h) 62 and peak post-mating gene expression in the spermatheca occurs at 24 h 23 . In the Ae. aegypti lower reproductive tract, more genes are differentially expressed at 24 h than at 6 h post-mating 41 .
To determine how Ae. aegypti spermathecae respond to mating and blood-feeding, we examined transcript abundance in the spermathecae at 6 and 24 h to uncover early changes. We also examined transcript abundance at 72 h, when oocytes for our mosquitoes had matured following a blood-meal and egg-laying/sperm release from storage commences.
We observed a large number of differences in transcript abundance patterns across the three time points, with a total of 797 differentially expressed transcripts in spermathecae from NBF females compared to 1,158 transcripts in the BF group ( The magnitude of the transcriptional response-the average log-fold change of the total number of differentially expressed transcripts compared to virgins-significantly differed between the NBF and BF groups in both up-and down-regulated transcripts across all timepoints (LM F = 8.89, P = 0.00015; Fig. 3E). In spermathecae from BF females, the magnitude of the transcriptional response was significantly higher at 24 h compared to 6 and 72 h for both up-regulated (P < 0.0001) and down-regulated transcripts (P = 0.001). However, the magnitude of the transcriptional response in spermathecae from NBF females was similar across the three evaluated timepoints for up-regulated transcripts (P > 0.05). Furthermore, a significant decrease in down-regulated transcripts occurred after the 6 h timepoint in NBF females (P = 0.01; Fig. 3E).
Blood-feeding after mating induces a transcriptional response in the spermathecae. Bloodfeeding by females is a requirement for anautogenous mosquito reproduction. However, the potential link between blood-feeding and the post-mating regulation of gene expression in reproductive tissues has not been explored. While mating induces a transcriptional response in reproductive tract tissues of female mosquitoes 23,31,41 , the extent to which this transcriptional response is mating-vs. blood-meal-induced is not yet known. The mating encounter site for Ae. aegypti is around the host-males intercept females as they come to blood feed 70 . Therefore, to examine post-mating changes in gene expression in a biologically-relevant situation, we identified spermathecal transcripts in females that had first mated and blood-fed shortly afterward (see Methods). As mating alone induces a transcriptional response in the spermathecae (Figs. 2A, 3A), our data do not permit determination of spermathecal transcripts regulated only by blood ingestion. Recurring blood-feeding events by mated females may regulate genes important for sperm viability or for egg-development, ovulation and/or egg-laying.
We compared transcript abundance between NBF and BF female spermatheca using: (1) the combined set of all timepoints for each female group, and (2) samples from each female type at each timepoint. We found only one differentially expressed gene when comparing the combined post-mating timepoints of BF spermatheca samples to NBF samples (FDR < 0.05; LogFC > 1; Figure S2C). This gene, phosphoenolpyruvate carboxy-kinase (PCK; AAEL000080), is associated with the generation of glucose from non-carbohydrate precursors, such as pyruvate and amino acids. The abundance pattern of this gene is significantly greater at 6 h in spermathecae from BF females and decreases at 24 and 72 h (Fig. 4A). Enzymes involved in glucose production, such as Scientific RepoRtS | (2020) 10:14899 | https://doi.org/10.1038/s41598-020-71904-z www.nature.com/scientificreports/ glucose-dehydrogenase in D. melanogaster, play a role in sperm nourishment and release from storage sites 71 . Thus, transcripts involved in catalytic reactions to produce energy in mosquito spermathecae might form a basis to discern the process of energy supply in the sperm storage organs. Next, we examined transcripts at each post-mating timepoint. We observed a decreasing number of transcripts that were up-regulated by blood-feeding from 6 h (21 transcripts) to 24 h (10 transcripts) (Fig. 4B, Figure S2J-L; File S5). By 72 h we did not identify differentially abundant transcripts when comparing mated BF and mated NBF females (Fig. 4B). At 6 h, transcripts with the highest change in expression were related to energy metabolism: enoyl-CoA hydratase (AAEL003993), D-amino-acid oxidase (AAEL000213), and phosphoserine aminotransferase (AAEL012578) (File S5). At 24 h, vitellogenin-A1 (AAEL010434), vitellogenic carboxypeptidase (AAEL006563), and the general odorant-binding protein 99a (AAEL005772) were the most significantly up-regulated transcripts (logFC > 3, FDR < 0.01). Although vitellogenesis is activated by blood-feeding, vitellogenin is up-regulated > 200 fold in the An. gambiae spermatheca post-mating 62 . Yolk proteins are also detected in the D. melanogaster sperm storage organs 32 . Many up-regulated transcripts in BF females are related to energy metabolism and amino acid biosynthesis, suggesting that blood-feeding might elicit production of molecules for sperm nourishment and energy provision. Previous studies of mated Ae. aegypti females have proposed that males transfer mRNA during mating-106 mRNAs have been identified as potentially male-derived 41,54 . We detected 6 of these transcripts in our differential expression analysis of spermathecal transcriptomes from NBF and BF females (File S6). Some mRNAs identified in the lower reproductive tract of mated females are derived from the male accessory glands 41,72 . Further, mammalian 73,74 and Drosophila 75 sperm contain a "pool" of RNAs. While sperm-specific RNAs have not been identified in Ae. aegypti, we wished to determine if the transcripts identified in our analysis might have originated from the testis or the accessory gland of the male reproductive tract. We compared the significantly up-regulated transcripts identified in the NBF and BF female groups to the male accessory gland and testis transcriptomes 72 . We found 11 transcripts that overlap with the testis transcriptome (0.56% of the total transcripts) and 36 that overlap with the accessory gland transcriptome (1.84%) (File S6). Of these genes, only one, the predicted headpeptide (AAEL024630), had no signal in the spermathecae of virgin females. However, we cannot rule out that the other 46 transcripts identified in the testes and/or accessory gland transcriptomes may have contributed to the signal we detected from the spermathecae of mated females. Further experimentation is required to determine the sex-specific expression of the testes/accessory gland transcripts also identified in our analysis.
Functional categories of differentially expressed transcripts. We examined the major functional categories represented in the differentially expressed spermathecal genes from NBF and BF females compared to virgin females for all combined time points. Significant differences in GO terms were observed in up-and down-regulated transcripts between the NBF and BF groups. In the up-regulated dataset, we found 48 significantly overrepresented GO terms (FDR < 0.05) in the NBF group (File S7), with the majority related to transport activity. Other terms are associated with ATP activity, pH regulation and peptidases (File S7). The most significant terms correspond to energy-coupled proton transmembrane transport against an electrochemical gradient, ATP hydrolysis-coupled proton transport, proton transmembrane transport, ATPase activity and active transmembrane transporter activity (Table 1). In the BF group, we found 24 overrepresented GO terms related to amino acid metabolic processes, carboxylic acid metabolism, inosine monophosphate (IMP) metabolism and biosynthesis, transport activity and peptidase activity (File S7). The most significant terms corresponded to ' de novo' IMP biosynthetic process, integral component of membrane, secondary active transmembrane transporter activity, IMP biosynthetic process and IMP metabolic process (Table 1).
In the down-regulated dataset, we found 14 and 58 significantly overrepresented GO terms (FDR < 0.05) for the NBF and BF groups, respectively. Enriched terms in the NBF group are related to membrane components, endoplasmic reticulum, peptidases, proteasome core complex, cell redox homeostasis and proteolysis (File S7). In the BF group, there was a significant down regulation of several functions (File S7), with the most significant terms corresponding to structural constituent of ribosome, ribosome, translation, ribosomal subunit and cytosolic large ribosomal subunit (Table 1).
We next performed GO enrichment analysis of differentially expressed transcripts from NBF and BF female groups compared to virgin females at each post-mating time-point. In the up-regulated dataset at 6 h, significant functions are only observed in the BF group (Table 2; File S8), with the most significant terms corresponding www.nature.com/scientificreports/ to IMP biosynthetic process, IMP metabolic process, oxoacid metabolic process and carboxylic acid metabolic process (Table 2). At 24 h, we found 83 and 26 significant terms for the NBF and BF female groups, respectively. Overrepresented terms in the NBF group are related to transport activity and ATP transport and metabolism (File S8). The most significant terms are energy-coupled proton transmembrane transport and ATP hydrolysis coupled proton transport (Table 2). In the BF group, several terms are related with transmembrane transporter activity and ion transport ( Table 2; File S8). Analysis of down-regulated transcripts showed statistically significant terms only at 24 and 72 h in both NBF and BF female groups (File S8). In NBF females, significant terms are related with endoplasmic reticulum and peptidases at 24 h (Table 2). An increase in the number of significant terms was observed at 72 h (File S8), with several related to proteolysis ( Table 2). In the BF group, the most significant terms at 24 h were RNA binding proton-transporting V-type ATPase, and translation activity ( Table 2). Significant GO terms increased at 72 h, related to ribosome, translation, protein metabolic process and structural molecule activity (Table 2). comparative expression patterns between the spermathecae and lower reproductive tract tissues. The Ae. aegypti lower reproductive tract comprises the bursa, oviduct, spermathecal vestibule, spermathecal ducts and the spermathecae 21,40 . To determine transcripts more abundant in the spermathecae compared to the other reproductive tract tissues, we generated an overview of transcripts identified in this study compared to transcripts identified in the lower reproductive tract by Alfonso-Parra et al. 41 . We developed a comparative analysis to identify differentially expressed spermathecal transcripts from NBF females with the differentially expressed transcripts identified in all lower reproductive tract tissues at 6 and 24 h post-mating.
A total of 116 differentially up-regulated transcripts were identified in both the NBF and BF datasets: 44 and 72 transcripts at 6 and 24 h, respectively ( Fig. 5; File S9). The most significant up-regulated transcripts correspond to tyrosine-protein kinase hopscotch (AAEL003619) of the JAK/STAT signaling pathway at 6 h, odorant binding protein OBP23 (AAEL006109), and brain chitinase and chia (AAEL002972) at 24 h ( Fig. 5; File S9). A total of 48 transcripts (18 at 6 h and 30 at 24 h) were found to be significantly down regulated in both studies ( Fig. 5; File S9). The most significant down-regulated transcripts correspond to adenylate cyclase type 8 (AAEL022948) at 6 h and farnesol dehydrogenase-like (AAEL009685) at 24 h. While this comparison identified several overlapping genes, the number of unique transcripts found in the former study 41 suggests that individual tissues of the Ae. aegypti female reproductive tract respond to mating. As mating induces changes in gene expression in reproductive tract tissues in addition to the spermathecae in other insects 69 , further analyses of Ae. aegypti female reproductive tract might uncover tissue-specific genes that are important in reproduction of this disease vector. comparison of the early and late mated female spermathecae transcriptome. A recent study by Pascini et al. 31 reported the spermathecae transcriptome in mated females at 7 d post-eclosion. We compared www.nature.com/scientificreports/ the differentially expressed genes in our datasets with those from the previous study (File S10), as these transcripts might represent long-term transcriptional changes in the spermathecae in response to mating. Seven days post-eclosion corresponds to ~ 6 d post-mating, as females require 24 h to reach sexual maturity and become receptive to mating 22 . Comparing the overall transcriptional response of the spermathecae in NBF females at 6, 24 and 72 h with that at 7 d post-eclosion, we found higher differentially expressed transcripts at each timepoint compared to the 7 d post-eclosion timepoint from Pascini et al. 31 . We also observed an overall decrease in gene expression levels in the NBF female group from 6 to 72 h, suggesting that the spermathecal transcriptional www.nature.com/scientificreports/ response decreases over time after mating. Moreover, 17 transcripts up-regulated at 6, 24 and 72 h were observed to be down-regulated at 7 d post-eclosion (File S10). Four up-regulated transcripts observed at 6 h in the NBF female group were found to be up-regulated at 7 d post-eclosion (File S10), suggesting that these genes return to high expression levels in the spermathecae nearly a week after mating. Pascini et al. 31 performed functional analysis on 8 genes found to be differentially regulated at 7 d posteclosion. Of these 8 genes, we found only one in our dataset: the putative potassium-dependent sodium-calcium exchanger AAEL004814. In NBF female spermathecae, AAEL004814 transcripts are significantly up-regulated at 72 h post-mating. However, AAEL004814 had a greater response to blood-feeding: AAEL004814 transcripts were significantly down-regulated at 6 h, and significantly up-regulated at 24 and 72 h in the BF female group (File S4). The regulation of AAEL004814 by blood-ingestion is consistent with its role in oocyte development 31 . conclusions Immediately after mating ends, female insects undergo physiological and behavioral changes necessary for the production of progeny 15,16 . While several post-mating changes in female insects are stimulated after receipt of male seminal proteins 15 , spermathecal products are also required for optimal fertility 24,27,35 , suggesting that male and female-derived molecules might coordinate molecular activities to maximize reproductive success. Postmating modification of female reproductive tract tissue 61 or sperm in storage 72 occur within hours of mating, and some seminal proteins that mediate sperm storage or ovulation require female contributions for optimal activity 76 or undergo proteolytic cleavage in the sperm storage organs that are necessary to elicit their effects 26,77 . Mating responsive genes may be required for post-mating tissue or sperm modifications to occur and/or to mediate the effects of seminal proteins within the female reproductive tract.
The identification of reproductive molecules that are required for fertility is an important endeavor in developing desperately needed new strategies for vector control. Elucidating changes in mRNA expression after mating and blood-feeding in female Ae. aegypti spermathecae can identify genes that are important for sperm storage, and that thus can serve to monitor female mating physiology or provide novel molecular targets to control this important vector in nature. Unlike previous transcriptome studies in Ae. aegypti reproductive tissues that identified differentially expressed genes in mated, sugar-fed females, this study examined transcriptional change in www.nature.com/scientificreports/ response to both mating and blood-feeding in the spermathecae. We demonstrate that the spermathecal transcriptional response to mating affects regulation of transcripts related to osmotic balance, membrane transport and oxidoreduction. Furthermore, blood-feeding induces greater transcriptional change in the spermathecae compared to mating alone. Our results reveal the complexity of gene regulation in the Ae. aegypti spermathecae in response to mating and a blood-meal in the immediate days after mating when females undergo physiological and behavioral changes that are necessary for the production of progeny.

Methods
Mosquitoes. Thai  Mating and spermathecae collection. We observed all matings by placing one female and three males into an 8 L bucket until a copulation occurred (defined as the engagement of male-female genitalia for ≥ 10 s 19,78 ). The pair were subsequently removed, and the male discarded. Females were grouped into 20 min mating intervals. A subset of females were blood-fed on the arm of a volunteer immediately following a mating interval; only fully engorged females were considered blood-fed. Subsequent to mating and blood-feeding, females were housed in 500 ml cups and maintained in the environmental chamber until spermathecal dissection. Bloodfeeding was approved by the Comité de Bioética Sede de Investigación Universitaria (Universidad de Antioquia); informed consent was obtained from all volunteers who participated in the study. All methods were performed in accordance with the relevant guidelines and regulations. Spermathecal samples (from 3 biological replicates) were taken from (1) virgin females, (2) mated females and (3) mated + blood-fed females at 6, 24, and 72 h post-mating (Fig. 1). Females were anesthetized on ice and spermathecae dissected in sterile, 1X phosphate-buffered saline. Fat bodies and other associated tissues were carefully removed. The sample corresponded to a pool of 25 spermathecae (three complete reservoirs and ducts). Dissected spermathecae were placed into 100 μl QIAzol reagent (Qiagen, Hilden, Germany), centrifuged at maximum speed for 1 min and stored at − 80 °C until mRNA extraction.
RnA extraction and library preparation. RNA extraction was based on a chloroform/isopropanol isolation and ethanol precipitation protocol following manufacturer's instructions. GlycoBlue (Thermo Fisher Scientific, Waltham, USA) was added as a co-precipitant prior to ethanol precipitation. RNA concentration was measured using a Qubit spectrophotometer (Invitrogen, Grand Island, NY). At least 10 ng total RNA was used as input for polyA enrichment of mRNA. cDNA libraries were prepared using the NEBNext Ultra II Directional RNA Library Preparation kit for Illumina sequencing following the protocol for polyA enrichment with the NEBNext Poly(A) mRNA Magnetic Isolation kit (New England Biolabs, Ipswich, MA). Library preparation was performed at the Cornell University Transcriptional Regulation and Expression facility (Ithaca, NY).
Sequencing and data analysis. Samples were sequenced using the NextSeq500 platform with 85nt single end reads at the Cornell Biotechnology Resource Center (Ithaca, NY). A minimum 20 M raw reads were generated per sample. Reads were processed using trim-galore (Babraham Institute) to filter low quality reads, filter noisy short fragment, and filter adapter sequences. Filtered reads were then aligned to the Ae. aegypti NCBI genome release (AaegL5.0) with RefSeq annotations using STAR 79 . A matrix of raw read counts for each annotated genomic feature was extracted. Genes with low read counts were filtered from the matrix (CPM ≤ 5), and RUVseq 80 was used to identify sources of erroneous variation using residuals (k = 3), which were included in a design matrix that was entered into edgeR 81 for differential abundance testing. Filtered reads were then aligned to the Ae. aegypti NCBI genome release (AaegL5.0) with RefSeq annotations using STAR. A matrix of raw read counts for each annotated genomic feature was extracted. Upon sample QC, we found that one replicate for virgin and one for the 24 h time point in BF females contained low quality data and were thus removed from downstream analysis. A total of 8,560 transcripts passed our quality control filters: 8,223 annotated and 337 unannotated transcripts on the VectorBase bioinformatics database (https ://www.vecto rbase .org). To perform gene-wise differential testing, we used a quasi-likelihood model and fit generalized linear models, and ultimately performed F-test to identify differentially abundant genes. To perform Gene Ontology (GO) enrichment analysis, we used ontology terms that were generated by Degner et al. 15 , and performed enrichment tests using the GOseq 82 package in R. The R analysis script is available with this manuscript as Supplementary File S11.

Data availability
Data are available on NCBI under Accession Number PRJNA612326.