Effects of temperature on transcriptome and cuticular hydrocarbon expression in ecologically differentiated populations of desert Drosophila

Abstract We assessed the effects of temperature differences on gene expression using whole‐transcriptome microarrays and cuticular hydrocarbon variation in populations of cactophilic Drosophila mojavensis. Four populations from Baja California and mainland Mexico and Arizona were each reared on two different host cacti, reared to sexual maturity on laboratory media, and adults were exposed for 12 hr to 15, 25, or 35°C. Temperature differences influenced the expression of 3,294 genes, while population differences and host plants affected >2,400 each in adult flies. Enriched, functionally related groups of genes whose expression changed at high temperatures included heat response genes, as well as genes affecting chromatin structure. Gene expression differences between mainland and peninsular populations included genes involved in metabolism of secondary compounds, mitochondrial activity, and tRNA synthases. Flies reared on the ancestral host plant, pitaya agria cactus, showed upregulation of genes involved in metabolism, while flies reared on organ pipe cactus had higher expression of DNA repair and chromatin remodeling genes. Population × environment (G × E) interactions had widespread effects on the transcriptome where population × temperature interactions affected the expression of >5,000 orthologs, and there were >4,000 orthologs that showed temperature × host plant interactions. Adults exposed to 35°C had lower amounts of most cuticular hydrocarbons than those exposed to 15 or 25°C, including abundant unsaturated alkadienes. For insects adapted to different host plants and climatic regimes, our results suggest that temperature shifts associated with climate change have large and significant effects on transcriptomes of genetically differentiated natural populations.


| INTRODUCTION
Understanding current responses to increasing temperatures in natural populations can help to predict future responses to climate change (Hoffmann & Blows, 1994;Williams, Elizabeth, & Samantha, 2003) because temperature is the primary environmental factor determining most organisms' distribution and abundance (Cossins & Bowler, 1987).
Warm deserts and the communities that inhabit them are particularly well suited to a molecule-to-ecosystem approach to understanding mechanisms of temperature adaptation because these extreme conditions directly impact organism survival and reproduction, which in turn are directly influenced by complex suites of behavioral, physiological, and biochemical traits. The arid lands of the southwestern USA and northwestern Mexico are predicted to become warmer and drier in the future due to climatic shifts (Cook & Seager, 2013). Thus, it is crucial to understand how organisms are likely to respond to increased desertification associated with climate change, and an important first step is to understand how current biological communities function in warm, arid environments.
Although many studies have examined genetic variation, physiology, and ecology of desert organisms, our understanding of genomelevel adaptation to desert regions is in its early stages (Matzkin, 2012;Matzkin & Markow, 2009;Rajpurohit et al., 2013). The first desert organism with a sequenced genome was Drosophila mojavensis (Gilbert, 2007), and the genomes of its close relatives, D. arizonae (Lohse, Clarke, Ritchie, & Etges, 2015) and D. buzzatii (Guillén et al., 2015), are now available for comparative study. These genomes provide powerful models for an integrated understanding of resistance to the harsh temperature conditions of arid lands through analysis of population genomic divergence and whole-transcriptome responses to extreme desert conditions. Here, we focus on D. mojavensis, a member of the well-studied cactus-yeast-Drosophila model system, with a rich ecological background for genomic studies (Barker & Starmer, 1982;Barker, Starmer, & MacIntyre, 1990) because it has undergone adaptation to different host plants (Etges, 1990(Etges, , 1993Etges, de Oliveira, Noor, & Ritchie, 2010), and diverged into a variety of ecologically distinct desert habitats (Etges, Johnson, Duncan, Huckins, & Heed, 1999;Heed, 1982;Heed & Mangan, 1986). Furthermore, geographically isolated populations in Baja California and mainland Mexico and Arizona separated by the Gulf of California are considered incipient species due to the evolution of behavioral differences caused by courtship song and cuticular hydrocarbon differences (Etges & Ahrens, 2001;Etges, Over, de Oliveira, & Ritchie, 2006;Wasserman & Koepfer, 1977). Thus, D. mojavensis is an excellent model for understanding genomic responses and adaptation to desert environments.
In the deserts and arid lands of the Americas, about half of the ca 100 species in the D. repleta group carry out their life cycles in fermenting cactus tissues (Heed, 1982;Wasserman, 1992). The three species, D. mojavensis and its two closest relatives, D. arizonae and D. navojoa, form a monophyletic group endemic to northwestern Mexico and the southwestern USA (Ruiz, Heed, & Wasserman, 1990). Across its geographical range, D. mojavensis uses several host cacti including pitaya agria cactus, Stenocereus gummosus, in Baja California and organ pipe, S. thurberi, and cina cactus, S. alamosensis in mainland Mexico and Arizona (Etges et al., 1999;Heed & Mangan, 1986). Isolation from mainland D. arizonae occurred ca 1.3 mya after tectonic drift of present-day peninsular Baja California from southwestern Mexico (Gastil, Phillips, & Allison, 1975), and gene exchange between D. arizonae and D. mojavensis continued until ca 250 kya based on genomic divergence estimates and ages of chromosomal inversions (Lohse et al., 2015). This was about the same time that D. mojavensis colonized mainland Mexico from Baja California (Smith, Lohse, Etges, & Ritchie, 2012). Other host cacti are used by D. mojavensis in southern California, but in this study we focused on the well-studied Baja California and mainland Mexico populations.
A previous study of the effects of desiccation exposure on adult D. mojavensis showed striking effects of host cactus, age, and population differences on survivorship, water loss, cuticular hydrocarbon variation, and transcriptome variation (Rajpurohit et al., 2013). Taken together, we hypothesized that transcriptomic and cuticular hydrocarbon responses to temperature variation should also vary among populations of D. mojavensis, and be influenced by exposure to different host cactus substrates during pre-adult stages. Therefore, we assessed whole-transcriptome responses to temperature variation in multiple populations from different biotic regions of the Sonoran Desert reared on different host plants. We assessed responses in cuticular hydrocarbon composition to temperature differences because of their roles in mate choice (Etges & Ahrens, 2001) and presumed barriers to transcuticular water loss (Gibbs & Rajpurohit, 2010).

| Population collections and husbandry
We collected D. mojavensis from four locations, two on each side of the Gulf of California: San Quintin and Punta Prieta from Baja California, and Punta Onah, Sonora and Organ Pipe National Monument, Arizona (Table 1; see Etges et al., 2010 for mapped locations). Flies were collected directly from cactus rots by aspiration and by sweep-netting over baits in the field, returned to the laboratory, and cultured on banana food (Brazner & Etges, 1993) in 35-mL shell vials at room temperature (20-22°C).
After several generations in the laboratory, thousands of adult flies from each population were introduced into 12,720 cm 3 Plexiglas ® population cages and allowed to mate for 7 to 10 days. For each of the four populations, eggs were collected from the cages and distributed into six 250-ml milk bottles containing banana food. These cultures were started at medium larval densities to minimize nutritionally caused maternal/environmental effects caused by mass rearing in vials and maintained in an incubator that cycled from 27 to 17°C on a 14-hr:10-hr LD photoperiod. All eclosed bottle-reared adults were separated by sex daily and aged to maturity on laboratory food in vials in the incubator.
Sexually mature adults of each sex were counted into groups of 100 and then introduced together into population-specific oviposition chambers. Adults were allowed to mate for several days, and fresh agar-cactus oviposition media was used to collect eggs for 10 hr in 5.5-cm-wide Petri dishes. Eggs from each chamber were washed in deionized water, 70% ethanol, sterile deionized water, counted into groups of 200, transferred to a 1-cm 2 piece of sterilized filter paper, and placed on pieces of fermenting cactus in the same incubator as described above. Fermenting cactus cultures were initiated in plugged 250-ml bottles with a 5.5-cm-diameter piece of filter paper covering 75 g of aquarium gravel. Bottles were autoclaved and cooled, 75 g of either agria or organ pipe tissues were placed in the bottles, and the bottles were autoclaved again for 8 min at low pressure. After cooling to room temperature, each culture was inoculated separately by injection with a hypodermic needle with 0.5 ml of an aqueous mixture of a pectolytic bacterium, Erwinia cacticida (Alcorn et al., 1991), and 1.0 ml of an aqueous solution of seven yeast species found in natural cactus rots (Starmer, 1982): Dipodascus starmeri, Candida sonorensis, Candida valida, Starmera amethionina, Pichia cactophila, Pichia mexicana, and Sporopachydermia cereana. The inoculant solutions were prepared with sterile water and a full inoculating loop of each species of bacteria or yeasts. After 1-2 days, all unhatched eggs were counted in order to calculate egg to adult viability and development time. All eclosed adults from each replicate culture bottle were counted daily under CO 2 anesthesia, separated by sex, and kept on banana food in vials.
All cactus culture bottles and vials containing adults were housed in an incubator at 25°C with a 14-hr:10-hr LD cycle until the experiments began.

| Temperature treatments
Groups of 30 male or female adults were placed into separate 35ml glass vials containing banana food and exposed to 15, 25, or 35°C for 12 hr. This time period was chosen to reflect half of a diurnal temperature cycle, and the range of temperatures were chosen to approximate average low and high temperatures in nature (Gibbs, Perkins, & Markow, 2003). After exposure, females were flash-frozen in liquid nitrogen and stored at −80°C prior to RNA extraction. The experiments always started at 8:00 a.m., so flies were always exposed to the day/light cycle in the incubators. Only 7-day-old females were used for RNA extraction and microarray analysis because of the large number of samples, but both females and males exposed to the temperature treatments were included in the CHC analysis. A total of 24 treatments (four populations × two cacti × three temperatures) were included with four replications resulting in 96 groups of females for gene expression assays. The experiment began in February 2009, ca 12-18 months after the flies were captured in the field.

| Analysis of epicuticular hydrocarbon variation (CHCs)
After exposure to each of the temperature treatments, 8-day-old female and male adults were frozen and stored at −20°C until CHC extraction. Total CHCs were extracted by rinsing individual adults in hexane for 20 min in a 300μl glass vial insert (Microliter Analytical Supplies, Suwanee, GA), evaporating off all solvent in a 40°C heating block, and then freezing each sample at −20°C. Individual CHC extracts were redissolved in 5 μl of heptane containing a known amount of docosane (C 22 ) as an internal standard. Each sample was subjected to capillary gas-liquid chromatography using an automated Shimadzu GC-17A (Shimadzu Scientific Instruments, Columbia, MD) fitted with a 15-m (ID = 0.22 mm) Rtx-5 fused-silica column (Restek Corporation, Bellefonte, PA). Injector and detector temperatures were set at 290 and 345°C, respectively, with the injector port in split mode (3:1 ratio), and the column was heated from 200 to 345°C at 15°C/min holding at 345°C for 4 min.
CHC amounts were quantified in all flies by quantifying peak integrations using Class VP 4.2 software provided by Shimadzu with C 22 amounts used as an internal standard and expressed as nanograms/ fly (Etges & Ahrens, 2001;Etges & Jackson, 2001;Stennett & Etges, 1997). All CHC data were log 10 -transformed to improve normality.
T A B L E 1 Origins of the four populations of Drosophila mojavensis in this study, the biotic regions of the desert they are located in (Etges et al., 1999), and numbers of flies used to establish laboratory populations. All flies were collected over banana baits in nature unless otherwise noted. We sampled five adults for each combination of population, sex, cactus, and temperature, and MANOVA was used to estimate significant sources of variation using the model: where μ is the grand mean, P i is the effect of population (two Baja California vs. two mainland populations), S j is the effect sex, C k is the effect of host cactus, T l is the effect of temperature, I P×S is the interaction between population and sex, I P×C is the interaction between population and cactus, I P×T is the interaction between population and temperature, I S×C is the interaction between sex and cactus, I C×T is the interaction between cactus and temperature, I P×S×C is the in-

| Statistical analysis
Our experimental design consisted of 96 hybridizations that were analyzed with a mixed-model ANOVA, PROC MIXED (SAS-Institute,

2004b):
where μ is the grand mean, R i is the effect of geographic region (Baja California vs. the mainland), P j is the effect of population nested within regions, C j is the effect of host cactus, T k is the effect of temperature, I R×C is the interaction between region and cactus, I R×T is the interaction between region and temperature treatment, I P×C is the interaction between population and cactus nested within region, and E ijkl is the error term. To correct for multiple comparisons, we calculated false discovery rates, FDR, for the overall ANOVA data and for all pairwise comparisons between treatments (Benjamini & Hochberg, 1995). All comparisons with FDR < 0.01 were included in our analyses.

| Bioinformatic analysis
We used DAVID 6.7 and 6.8 beta (Huang, Sherman, & Lempicki, 2009) to identify functional categories that were enriched as a function of each main experimental variable. We used DAVID's default setting of medium stringency to minimize exclusion of potentially interesting gene ontology (GO) terms. Clusters with enrichment scores <1.3, corresponding to p > .05, were excluded from further analyses. We submitted GO terms to GO module (Yang, Li, Lee, & Lussier, 2011) to identify false positives included because of the hierarchical nature of the GO terms. After identifying differentially expressed genes and enriched gene clusters, we used FlyMine (Lyne et al., 2007) to identify orthologous genes in D. melanogaster. FlyMine and FlyBase (Attrill et al., 2016) were used to explore gene functions and tissue-specific expression patterns.

| RESULTS
A total of 7213 adult D. mojavensis were reared on both host cacti (n = 48 cultures). Nested ANOVAs of both egg to adult development time (DEVT) and viability (Table S1) revealed similar life history patterns to those in previous studies (Etges, 1990;Etges et al., 2010).  Figure   S1). All four populations showed increases in DEVT when reared on agria vs. organ pipe cactus, due most likely to the more rapid tissue breakdown of agria tissues at 25°C-in all previous studies, cultures were grown on a 14-hr:10-hr LD photoperiod and 27:17°C temperature cycle, a mean daily temperature of 22.8°C.
Baja California populations expressed higher egg to adult viability than mainland populations, 82.3 vs. 75.1, two-tailed t-test, p = .016, and there were significant differences between populations; F = 3.37, p = .0276, and cactus substrates; F = 9.23, p = .0042, but these were not apparent in the nested model (Table S1). Agria cactus decreased viability in both mainland populations, but both Baja California populations showed increased viability on both cacti ( Figure S1) similar to previous studies (Etges et al., 2010).

| Epicuticular hydrocarbon variation
MANOVA revealed that every main and interaction effect in the complete model was significant, p < .0001 (Table 2) revealing the sensitivity of CHC expression due to population differentiation, gender, rearing substrates, temperature, and their interactions. We did not employ a nested model to simplify analysis because large-scale regional differences in CHCs have already been documented (Etges & Ahrens, 2001;Etges & Tripodi, 2008). Principal components analysis revealed eight PCs that accounted for 86.1 percent of the variation in the data (Table S2). ANOVA was performed using the same full factorial model to reveal the nature of each of the eight PCs (Table 3).
Population, cactus, and sex influenced each of the first three PCs, and temperature influenced six of eight PCs. Each of the eight PCs showed significant model effect interactions with temperature revealing a pervasive influence of these treatment effects on different covarying combinations of CHCs in this study.
Post hoc multiple comparisons of the effects of temperature on each of the 31 CHC components and total CHCs per fly revealed that almost all CHCs were reduced in abundance at 35°C. A common pattern for many CHCs was for flies exposed to 15 and 25°C, CHC amounts were not significantly different, but greater than amounts observed at 35°C indicating decreased CHC production from 25 to 35°C (Table 4), consistent with some earlier observations (Markow & Toolson, 1990). Further, some CHC amounts at 15°C were slightly greater than at 25°C. Several CHC components including 7-and 9-hentricontene, 7,27-pentatricontadiene, C36a, C36b, 9,27-heptatricontadiene, C38, and C40 were higher in abundance at 25 °C, as were total CHCs per fly. As most of these CHCs are monoenes or alkadienes not expected to decrease cuticular permeability (Gibbs, 1998(Gibbs, , 2002, it was unexpected to observe decreased abundances of these CHCs at 35°C in just 12 hr of exposure. The most abundant alkane, 2-methyloctacosane, should be best at increasing desiccation resistance because it is a saturated hydrocarbon, or harder wax, but it also decreased in abundance at 35°C. Thus, there was little or no evidence that responses in observed CHC production at 35°C were associated with altered transcuticular water loss. Of those CHCs previously associated with differences in sexual attractiveness (Etges & Tripodi, 2008), both C 34 dienes shared this pattern of lowered expression at 35°C, and only one of the C 37 dienes, 9,27-heptatricontadiene, showed a significant response to temperature where higher amounts were expressed at 25°C (Table 4). This contrasted with increased total C 37 abundance (35-methylhexatricont-10-ene, 9,27-heptatricontadiene, 8,28-heptatricontadiene, and 14-, 16-, and 12-hexatricontene) at 31 or 34°C in laboratory food-reared flies (Markow & Toolson, 1990;Reidy, Toolson, & Markow, 1991 Eigenvalues and percent of total variance for each PC are provided in Table S2. but significant effect of cactus on CHC variation. Assessing host cactus effects at each temperature required splitting the data into six groups revealing how much of the variation in CV2 was due to organ pipe-reared adults at 35°C (Figure 1). All of the temperature-cactus centroids were significantly different from each other (p < .001) except for organ pipe-reared flies at 15 and 25°C. Thus, effects of temperature variation on 8-day-old D. mojavensis CHCs depended on the host cacti they were reared on from egg to eclosion even though they were held on laboratory food after eclosion.

| Patterns of gene expression
Each of the main effects affected the transcription of >1,400 genes, with temperature affecting 3,294 genes. Overall, nearly 80% of all predicted genes were differentially expressed in response to variation in temperature, region, population, host plant, or interactions between these (Table 5). Over 3,000 genes were differentially expressed, mostly as a result of pairwise and three-way interactions between main effects. Because DAVID is limited to sets of <3,000 genes, we were forced to restrict our analyses to pairwise comparisons for interaction effects where differences in expression were >1.5-fold in most cases.

| Temperature
Temperature affected the expression of 3,294 genes, including 2,490 with some annotation, that is, D. melanogaster orthologs. When all such genes were analyzed using DAVID, we identified five clusters of significantly overrepresented gene ontology (GO) terms ( P-450 proteins. The latter group included orthologs of four of the six "Halloween" genes (Rewitz, O'Connor, & Gilbert, 2007) in the biosynthetic pathway for 20-hydroxyecdysone, that is, disembodied, shadow, shade, and phantom. All were upregulated at 15 vs. 25 or 35°C (Table   S4). Expression of desat-2 also decreased as temperature increased, as expected for maintenance of proper membrane fluidity. This may also explain why lower amounts of CHCs were observed at higher temperatures, most of which have double bonds produced by desaturation of hydrocarbon precursors, for example, myristic acid (C14:0) by desat-2 in D. melanogaster (Yew & Chung, 2015). The gene clustering term "secondary metabolites biosynthesis, transport, and catabolism" was also included in this cluster, indicating temperature-related shifts in metabolism of host plant compounds. The second cluster was annotated as peptidase inhibitors, with 10 serpin genes of a total of 32 genes. The third and fourth clusters included genes associated with carbohydrate, amino acid, and lipid metabolism. The final cluster was associated with transcription and translation, including many ribosomal proteins and several genes in the His3 and His4 families (Table 6).
To explore transcriptome changes in more detail, we examined progressive changes in gene expression as temperature increased. A total of 894 genes differed in gene expression from groups of adults reared at 15 to 25°C, where 254 genes increased in expression and 640 genes were downregulated ( with endoplasmic reticulum indicating increased intracellular protein processing and storage. For downregulated genes, overrepresented GO categories included clusters for serpins, mannosidases, and proteoglycan recognition proteins ( Table 6). Members of these genes have been implicated in immune function, suggesting a reduced immune response at 25°C. We also examined the genes with the greatest changes in expression between 15 and 25°C and found that genes for several vitelline membrane proteins increased in expression, while several chorion protein genes decreased in expression (Table S4). In addition, there were 12 candidate T A B L E 6 Gene ontology and enrichment for the effects of temperature in this study. All functional clustering was based on genes with FDR p < .01 and > 1.5-fold change for each treatment effect and interaction To further explore temperature effects on gene expression, we examined those genes with the greatest overall change in expression. The 25 genes with the greatest increase in expression from 25 to 35°C included seven heat-shock proteins, including six of the 12 genes with the greatest increase in expression (Table S4) with chorion formation, was also in this group. Together with decreased expression of CG10407, decreases in lipid metabolism and egg production signaled reduced reproductive effort when adults were exposed to 35°C.
We also examined genes that showed significant shifts in expression (p < .01, fold change >1.5×) with comparisons involving all three temperature treatments. The comparisons (

| Regional differences between mainland Sonora/Arizona and Baja California
While we included just two populations from each region, we found 1,432 genes were differentially expressed in adult female D. mojavensis from Sonora/Arizona vs. Baja California populations (Table 5).
Approximately equal numbers of genes were upregulated and downregulated in these two regions ( Cyt-P450 genes, while the second reflected differences in transcriptional activity, including aminoacyl-tRNA synthetases. The other clusters were enriched in endopeptidases, including serpin genes, extracellular recognition proteins including peptidoglycan recognition proteins and peritrophins, and five small heat-shock proteins from the Hsp20 gene family (Table S6).
A total of 745 genes were expressed at higher levels in flies from Baja California than Sonora and Arizona of which 482 had D. melanogaster orthologs (Table 7). DAVID (Huang et al., 2009) revealed three clusters, including one enriched in GO categories related to peptidase activity, including four Jonah family genes. The Jonah genes are most highly expressed in the digestive tract, but other genes in this cluster are expressed in a variety of larval and adult tissues. The second cluster included several electron transport chain genes, as well as phm and 15 other Cyt-P450 genes. Interestingly, these were different P450 genes than those that were upregulated in mainland populations (above) except for genes orthologous to Cyp6a9 and Cyp6a21 that corresponded to more than one Dmoj ID suggesting problems with annotation (Table S6). The third cluster contained 10 peptidase genes, six of which are highly expressed in the digestive tract, while three have their greatest expression in the adult head (Attrill et al., 2016).

| Population differences
For the two populations each from Sonora/Arizona and Baja California, with populations nested within regions in our statistical model, a total of 2,433 genes were differentially expressed in pairwise comparisons between populations. This was reduced to 745 genes with FDR p < .01 and >1.5× fold change (Table 5). Not surprisingly, the majority of significant pairwise differences (84%) involved populations from different regions. Within regions, 587 genes differed in expression between the two mainland populations vs. 153 genes (83 annotated) between the two Baja California populations (Table 7). Cluster analyses using DAVID yielded results similar to those for region (Table S7), where the five most significant clusters (of 10) were generally categorized as peptidase activity, secondary metabolism, tRNA processing, peptidase inhibitor activity, and tRNA aminoacylation.
Comparisons of populations within each region, however, revealed significant differences in expression for different up-and downregulated gene clusters (  (Table 7).
T A B L E 7 Gene ontology and enrichment for the effects of region, populations, host cacti on gene expression differences in female Drosophila mojavensis in this study. All functional clustering was based on genes with FDR p < .01 for each treatment effect

Comparison No. genes (No. annotated) GO term Enrich score
Greater numbers of genes showed expression differences between the two mainland populations than between the two Baja California populations (  (Table 7). Thus, closer inspection of local population differences in gene expression revealed functional differences not seen in the broader regional comparisons.

| Host plant differences
Pre-adult rearing on different host cacti affected the expression of 2,457 genes of which 1,816 were annotated (FDR p < .01), but only six of these pairwise differences exceeded a fold change of 1.5× (Table 5), and all were overexpressed on agria vs. organ pipe cactus. Of these, five were annotated, with Dmoj_GI17776 and Dmoj_GI20850 involved with peptidase activity and Dmoj_GI11137 associated with transmembrane function in the major facilitator superfamily (MFS). This gene family includes single-polypeptide secondary carriers involved with transporting small solutes in response to chemiosmotic ion gradients (Mitchell et al., 2014). Dmoj_GI21891 and Dmoj_GI13378 were also overexpressed in agria-reared flies where the former has phospholipase activity and the latter is orthologous to flightin associated with adult muscle development in D. melanogaster (Table S8). Thus, for genes exceeding 1.5× fold changes, fermenting agria tissues increased expression of a small number of genes associated with ion membrane transport, protein and lipid catabolism, and muscle fiber development in adult D. mojavensis.
Because so many genes differed in expression due to cactus rearing substrates at FDR p < .01, but did not reach our 1.5× fold change threshold, we also investigated patterns of functional enrichment without this latter filter. Fifteen overrepresented clusters of GO terms were identified from genes there were differentially expressed at FDR p < .01 (Table S8), where the first cluster contained 137 peptidaserelated genes, and second included 41 Cyt-P450s and 43 other genes related to secondary metabolism including desat-2. Six clusters were associated with mitochondrial functions, particularly the electron transport chain. Three clusters were generally associated with extracellular processes, such as cellulose metabolism, chitinase, and extracellular matrix genes.
Striking differences were seen when genes were separated by their up-or downregulation on agria vs. organ pipe cactus. A total of 2,094 genes were expressed at higher levels on agria, while only 363 were upregulated on organ pipe cactus. Because genes upregulated on agria dominated the total pool, our clustering analysis was similar to that for all cactus substrate-affected genes. Eighteen clusters of GO terms were identified, including peptidases, secondary metabolism, and six mitochondrial function clusters (Table 7). In contrast, genes upregulated on organ pipe cactus were dominated by DNA repair (Table S8), replication, and chromatin structure, with five of six clusters relating to DNA functions. Twelve histone genes were upregulated in flies reared on organ pipe cactus. The remaining cluster of 39 ATP and purine nucleoside-binding proteins included two Hsp70 genes. The significance of these cactus-related gene expression differences is striking because they were caused by pre-adult rearing substrates; all RNA samples were extracted from 8-day-old adults reared on laboratory media, further evidence of carryover effects of pre-adult conditions into adulthood (Etges, de Oliveira, Rajpurohit, & Gibbs, 2016).

| Interactions between main effects
Clusters of differentially expressed genes identified by DAVID revealed detailed insights into the five pairwise interactions (Table 5) between regional, host plant, and temperature effects. Three-way interactions resulted in too many genes (> 3000) to be analyzed with DAVID, but for pairwise interaction terms including temperature, we concentrated on cactus and region interaction terms. There were 4,837 genes, 3,536 with some annotation, which differed in expression for Temperature × Host plant interactions resulting in eight enriched clusters (Table S9). The first included 95 proteases, especially Jonah family and trypsin genes. The second included 28 protease inhibitors, especially serpins. The third and fifth were enriched for
a Organ Pipe National Monument, Arizona and Punta Onah, Sonora, Mexico.
T A B L E 7 (Continued) Cyt-P450 and heat-shock genes, respectively. Each of these general clusters was found in one or both separate main effect analyses.
Inspection of functionally enriched gene clusters identified from up and downregulated genes for each significant contrast of cactus and temperature revealed a detailed portrait of temperature effects not seen in the overall analysis of Temperature × Host plant interactions (Table S9). Although these comparisons are not independent, they revealed consistent increases in expression of P-450, endopeptidase inhibitor, and juvenile hormone-binding protein gene clusters at low temperatures. At higher temperatures, Hsp20 and Hsp70 heatshock proteins, endopeptidase, and peptidoglycan recognition protein gene clusters were overexpressed, the latter involved with immune response and chitin binding (Table S10). A few enriched clusters were observed that were cactus specific such as α amylase-glycosyl hydrolase expression at 35°C and increased chorion gene expression in 15°C agria-reared flies.
For genes showing significant Region × Temperature interactions, DAVID recovered eight clusters of overrepresented GO terms. Several mirrored those found for temperature or region alone. The leading cluster represented a variety of metabolic processes, including secondary metabolism that involve metal ion binding, including three metallothionein genes, SOD3, desat-2, and carbohydrate metabolism genes (Table S11). Chromatin remodeling (ADD1) and Halloween (phm, dib) genes were also present. Other clusters were enriched for proteases, often expressed in the digestive system, and protease inhibitors. A small cluster (enrichment score 2.3) included five of the six heat-shock genes that were upregulated from 25 to 35°C. Upon closer inspection of all pairwise region-temperature differences, all of these Hsp20 family heat-shock genes were upregulated at 35°C in populations from both regions (Table S12).
For Region × Host cactus interactions, 3,999 unique genes showed significant differences in expression (FRD p < .01, > 1.5-fold ×), yet just 2711 were annotated (Tables 5, S13 and S14). The first cluster was enriched for 43 proteolysis/peptidase genes including trypsin encoding genes and the ortholog Ance, a member of the peptidase M2 family. There were also members of the M13 class of metalloproteases including neprilysins, involved with peptide hormone metabolism in Drosophila, and male and female reproduction including egg laying and sperm use in mated females (Sitnik et al., 2014). Some of these endopeptidase genes are expressed in the female reproductive tracts of D. mojavensis and D. arizonae (Kelleher & Markow, 2009;Kelleher & Pennington, 2009). Almost all were upregulated in Baja California populations and in flies reared on agria cactus (Figure 2). The second cluster included 47 genes enriched for P-450 activity and tetrapyrrole binding, and the third cluster was enriched for peptidoglycan function associated with microbial immunity (Tables S13 and S14).
Breaking the Region × Host cactus interactions into pairwise region-cactus comparisons revealed that upregulation of endopeptidases involved agria-reared flies from both regions (Table S13), and in one case Baja California flies reared on organ pipe cactus. In all cases, mainland populations showed upregulation of GO clusters enriched for secondary metabolites biosynthesis, transport, and catabolism, electron carrier activity, heme binding, or tetrapyrrole binding consistent with overall regional differences (Table 7). A number of other enriched clusters were noted that did not appear in the overall Region × Cactus term analyses were increased glutathione S-transferase (Gst) gene expression in mainland adults reared on both cacti, increased ATP production in the Baja California flies reared on agria vs. mainland flies reared on organ pipe cactus, and increased DNA repair and DNA synthesis in mainland adults reared on organ pipe cactus. Gst is involved with catabolism of toxic compounds in insects and shows nucleotide sequence and expression differences between populations of D. mojavensis and D. arizonae (Matzkin, 2008). Thus, the shift to organ pipe cactus in most mainland populations involved increases in secondary compound detoxification afforded by increased Gst and P-450 gene expression with a concordant increase in DNA repair systems.

| DISCUSSION
Temperature is thought to influence a majority of organismal bio- not included in this cluster, but GI24337, an ortholog of Droj2, is heatshock protein DnaJ included in gene cluster 2. In any case, heat-shock genes in D. mojavensis were upregulated well below lethal temperatures where many heat-shock proteins act as molecular chaperones to minimize thermal damage to cellular proteins, even though some damage is unavoidable. We also found evidence for increased expression of genes in ubiquitination pathways that would allow D. mojavensis to remove damaged proteins more rapidly at high temperatures (Zhou, Campbell, Stone, Mackay, & Anholt, 2012).
In addition to activation of heat-shock genes, an important part of the classical heat-shock response is the general downregulation of other genes. Our results also showed that chromatin remodeling affecting gene expression is an important component of the response to temperature. In genes whose expression increased from 25 to 35°C, the second significantly enriched cluster comprised genes associated with chromatin structure and modification. HDAC-6S, for example, upregulates genes associated with responses to abiotic stress (Cho, Griswold, Campbell, & Min, 2005). lid is an H3K4 methylase which is also required for H3 acetylation (Lloret-Llinares, Carré, Vaquero, de Olano, & Azorín, 2008), as is nej (Tie et al., 2009).
vig2 tends to increase H3K9me2 levels which are generally associated with a decrease in transcription (Gracheva, Dus, & Elgin, 2009).
Our transcriptome data also indicated changes in energy metabolism associated with temperature. From 25 to 35°C, genes involved in glycolysis and the pentose-phosphate shunt were downregulated.
Decreased expression of desat-2 was expected, as this would reduce membrane unsaturation and maintain membrane fluidity. Several other lipid biosynthesis genes were also highly downregulated (Table S4), while several lipases had greater expression. Together, these findings suggest that D. mojavensis shifts its fuel source from carbohydrates toward lipids as temperature increases. This shift may be affected by HDAC-6, which downregulates genes associated with glycolysis and biosynthesis (Cho et al., 2005). , but a few showed increased expression in mainland populations reared on agria (SO:AG). All were identified in an enriched functional gene cluster by DAVID (Huang et al., 2009) with an enrichment score = 12.6 (p < .0001; Table S11) 15 clusters enriched for GO terms included mitochondrial function, indicating that flies reared on agria cactus have higher metabolic rates.
DAVID also identified six clusters of overrepresented GO terms for genes upregulated on organ pipe cactus. The first included 25 genes associated with DNA repair, and all clusters were associated with aspects of nucleic acid metabolism (Table S8). Thus, flies reared on organ pipe cactus may be subjected to more DNA damage than those reared on agria.
A few broad categories of GO term clusters appeared repeatedly, especially in analyses of interaction effects. For example, clusters containing large numbers of proteases or protease inhibitors were found for two-way interactions between host plant and region or population (Table S13, Figure 2). The protease clusters often included trypsin or-  (Linder, Owers, & Promislow, 2008).
Flies from Baja California and Sonora/Arizona differed in the expression of numerous cytochrome P-450 genes, many of which are involved in detoxification and metabolism of secondary compounds in their host plants (Fogleman, Danielson, & MacIntyre, 1998;Frank & Fogleman, 1992). Agria and organ pipe cacti contain 25-45% triterpene glycosides and lipids to store carbohydrates and deter herbivory (Kircher, 1977(Kircher, , 1982

| Cuticular hydrocarbon expression
Short-term response to higher temperatures significantly decreased amounts of almost all CHC major components (Table 4) where all experimental model effects were highly significant (Table 2), and the downregulation of 12 CHC candidate desaturases and elongases.
Temperature differences resulted in overlapping CHC groups of cactus-reared flies, with only organ pipe-reared flies exposed to 35°C as adults displaying altered expression patterns (Figure 1)

| Climate change, desert ecosystems, and the responsive transcriptome
The Sonoran Desert and adjacent arid lands include well-studied biotic regions characterized by different plant and animal communities that have been shaped primarily by climatic differences (Shreve & Wiggins, 1964;Turner, Bowers, & Burgess, 1995;Turner & Brown, 1982 (Hastings, 1964;Hastings & Humphrey, 1969;Hastings & Turner, 1965), and all are within the species range of D. mojavensis.
Although the species assemblages in these biotic regions are derived from groups in the Miocene (8 mya), the Sonoran Desert developed ca 9 kya, with the current biotic regions evolving 4.5 kya (Betancourt, Van Devender, & Martin, 1990). Baja California separated from the Mexican mainland about 5.5 mya due to tectonic drift and has been in its present position for ca 1 mya (Murphy, 2006;Riddle, Hafner, Alexander, & Jaeger, 2000). Over this time frame, D. mojavensis populations expanded throughout these arid lands and now show significant genetic population structure associated with these biotic regions based on variation in chromosomal (Etges et al., 1999) and molecular data (Machado, Matzkin, Reed, & Markow, 2007;Smith et al., 2012).
Interpretations of differences in temperature adaptations between populations of D. mojavensis rest on these well-studied temperature-rainfall differences (Brown & Lowe, 1980). Although global climate models predict that the southwestern desert regions of North America will become warmer and drier in the next few decades (Alwelaie, Hoffman, & Saunders, 1995;Cook, Smerdon, Seager, & Coats, 2014), D. mojavensis may be buffered against the immediate impacts of higher temperatures as it has the highest thermal tolerance of any drosophilid studied to date. Larvae and adults survive shortterm exposure up to 44°C (Krebs, 1999;Krebs & Bettencourt, 1999;Stratman & Markow, 1998), and larvae do not display a significant heat-shock response until 38°C and above (Krebs & Bettencourt, 1999). Mainland adults continue to court, but not mate, up to 38°C (Patton & Krebs, 2001). This is especially relevant to understanding how D. mojavensis populations "over-summer" when temperatures are typically at their annual maxima, cactus rots are few, and fly populations are reduced to low levels (Rockwood-Sluss, Johnston, & Heed, 1973;R. H. Thomas, pers. comm.). Cactus microclimate temperatures in mainland Sonora and Arizona can exceed 42°C (Gibbs et al., 2003;Hastings & Humphrey, 1969), with average daily July temperatures >31°C in the Central Gulf Coast region of Baja California (Hastings, 1964). Increased expression of small heat-shock proteins in Sonoran/ Arizona adults vs. Baja California flies is consistent with cooler temperatures there because of the influence of Pacific coastal weather patterns. tRNA synthase expression was also higher in Sonora/ Arizona populations, suggesting higher overall transcription rates. In flies from Baja California, many electron transport chain and potential digestive proteases were upregulated, suggesting higher rates of metabolic turnover in flies from this region, consistent with findings of Rajpurohit et al. (2013).

| CONCLUSIONS
Uncovering genomewide patterns of transcriptional variation is only a first step in marshaling the data needed to predict the consequences of increasing temperature due to climate change in desert ecosystems. Proteomic and metabolomic changes, as well as whole organism fitness effects, will help to uncover the consequences of increasing temperatures to local adaptation and long-term population persistence. Both up-and downregulation of different gene networks influenced by population history and host plant effects emphasizes the need to place temperature-caused transcriptome variation in a broader ecological and biogeographical context. Expected upregulation of heat-shock proteins as well as chromatin remodeling genes associated with decreased transcription and reproductive effort showed the pervasive influences of temperature on stress responses, physiology, and fitness. While host plant use, temperature, and desiccation (Rajpurohit et al., 2013) had significant influences on CHC variation, seemingly nonadaptive decreases in CHC amounts and poorly understood temperature effects on sexual selection and isolation at higher temperatures will require more work to understand the multiple roles of CHCs in cuticle permeability and mate choice.
Understanding why agria cactus, the ancestral and preferred host, caused increases in expression of protein degradation and energy metabolism genes, including cyt-P450 genes, and organ pipe cactus increased expression of DNA repair and replication genes (Table 7) will also help to resolve the genetic basis of host plant adaptation in this oligophagic insect.

ACKNOWLEDGMENTS
We thank A. J. Marlon for laboratory wizardry, C. Ross for statistical assistance, K. C.