Wing morphology is related to host plants in cactophilic Drosophila gouveai and Drosophila antonietae ( Diptera , Drosophilidae )

A central issue in evolutionary biology is to understand the mechanisms promoting morphological evolution during speciation . In a previous study, we showed that the Neotropical cactophilic sibling species Drosophila gouveai and Drosophila antonietae can be reared in media prepared with their presumptive natural host plants ( Pilosocereus machrisis and Cereus hildmaniannus ) and that egg to adult viability is not independent of the cactus host. In the present study, we investigate the effects of ecological and genetic factors on interspeciﬁc divergence in wing morphology, in relation to the pattern of wing venation and phenotypic plasticity in D. gouveai and D. antonietae , by means of the comparative analysis of isofemale lines reared in the two cactus hosts. The species differed signiﬁcantly in wing size and shape, although speciﬁc differences were mainly localized in a particular portion of the wing. We detected signiﬁcant variation in form among lines, which was not independent of the breeding cactus, suggesting the presence of genetic variation for phenotypic plasticity and wing shape variation in both species. We discuss the results considering the plausible role of host plant use in the evolutionary history of cactophilic Drosophila inhabiting the arid zones of South America. Biological Journal of the Linnean Society , 2008, 95 , 655–665. evolution – phenotypic plasticity.


INTRODUCTION
A central issue in evolutionary biology is to understand the mechanisms promoting morphological evolution. Insect wings, probably causally related to the spectacular diversification of insects, are a particularly attractive system for comparative morphological studies (Carroll, Weatherbee & Langeland, 1995;Gibson, 1999). Specially, the wing of Drosophila is a suitable model for studies of morphological evolution in natural populations for several reasons (Bitner-Mathé & Klaczko, 1999a, b;Huey et al., 2000). First, many homologous landmarks can be identified and details of wing development are well understood (De Celis, 2003 and references therein). Second, it is a highly plastic organ in relation to a variety of environmental variables (David et al., 1994;Moreteau et al., 1998;Morin et al., 1999;Carreira et al., 2006). Third, wing morphology provides valuable phylogenetic information (Moraes et al., 2004a). Finally, different aspects of wing morphology are known to be the targets of natural selection Huey et al., 2000). Nevertheless, little is known about the evolutionary genetics of wing form. Some studies have found evidence consistent with optimizing selection acting on wing shape in natural populations (Weber, 1990;Gilchrist & Partridge, 2001), whereas others revealed high heritability for wing shape (Bitner-Mathé & Klaczko, 1999a, b;Moraes et al., 2004b). However, most of these studies have an important limitation; the ecology (breeding sites) of the species studied is poorly known. Species whose breeding sites are known have the added advantage of simulating appropriate ecological conditions in the laboratory as well as the possibility of extrapolating the results to natural conditions.
The study of phenotypic plasticity has become an important issue in ecological genetics (Conner & Hartl, 2004). Phenotypic plasticity refers to the property of a genotype to produce different phenotypes in response to an environmental change (Via et al., 1995;Schlichting & Pigliucci, 1998) and is an important phenomenon that has received increasing attention in studies of morphological variation (Debat et al., 2003). Regardless of its intrinsic relevance in evolution (Via et al., 1995;Schlichting & Pigliucci, 1998;Debat et al., 2003), the empirical study of the genetic basis of phenotypic plasticity has been difficult due to the polygenic nature of most traits studied and the limited knowledge of their usual environmental-dependent mechanisms of expression (Ungerer et al., 2003). In addition, different genotypes may vary in their response to a change in environmental conditions and, thus, reaction norms may take several forms (Gomulkiewicz & Kirkpatrick, 1992;Gavrilets & Scheiner, 1993). Such variation in reaction norms (or environmental sensibility) may be demonstrated in the form of a genotype-environment interaction (Falconer, 1990;Lynch & Walsh, 1998). It is well known that wing size and shape may respond to environmental variation in complex ways (Weber, 1990;Thomas & Barker, 1993;Bitner-Mathé & Klaczko, 1999b), suggesting that reaction norms may be part of an adaptive response.
The Drosophila buzzatii cluster (group repleta, subgroup mulleri, Wasserman, 1992) represents an ideal model for ecological and genetic studies because their strategies to cope with microenvironmental heterogeneity are mainly characterized by the use of well known cactus hosts (Hasson et al., 1992;Fanara, Fontdevila & Hasson, 1999). Recent studies in D. buzzatii (Patterson & Wheeler) and Drosophila koepferae (Fontdevila & Wasserman), the two basal members of the cluster (Manfrin & Sene, 2006), revealed changes in wing morphology that apparently were driven by the cactus hosts (Carreira et al., 2006). However, in the other five members, there are no studies linking wing morphology and host plant use (Manfrin & Sene, 2006).
Drosophila gouveai (Tidon-Sklorz & Sene) and Drosophila antonietae (Tidon-Sklorz & Sene) are two closely-related species (Manfrin, de Brito & Sene 2001;Manfrin & Sene, 2006) that also belong to the D. buzzatii cluster. The former is distributed from mid western Brazil to the Southern boundary of the Cerrado Domain (Tidon-Sklorz & Sene, 2001). In the southern section of its geographical range, D. gouveai is found in sandstone table hills associated with the columnar cactus Pilosocereus machrisis (Dawson) (Monteiro & Sene, 1995;Tidon-Sklorz & Sene, 2001;Moraes & Sene, 2002). In the north and north-east portions of its distribution, it occurs in large areas with great cactus diversity, including species of the genera Pilosocereus and Cereus (Pereira, Vilela & Sene, 1983). Drosophila antonietae (Tidon-Sklorz & Sene, 2001) is found in South and South-eastern Brazil and North of the Eastern border of the Argentinian Chaco, where it is associated with the columnar Cereus hildmaniannus (Schum) (Manfrin & Sene, 2006). Current evidence suggests that the southern limit of the distribution of D. gouveai is very close to the northern limit of D. antonietae. In this region, C. hildmaniannus and P. machrisis grow close together inside an east-west corridor of approximately 200 km (Manfrin & Sene, 2006).
In the present study, we investigated the effect that growing in different host cacti may have on wing morphology in D. gouveai and D. antonietae. We also identified the main sources of phenotypic variation by analysing wing morphology in isofemale lines reared in media prepared with the putative cactus hosts of each species.

MATERIAL AND METHODS
The present study is based on the analysis of isofemale lines founded with the progeny of wild inseminated females collected in natural populations of D. gouveai and D. antonietae (for details of the sites of collection, see Soto et al., 2007). Flies were reared in identical conditions for three generations in bottles with 30 mL of standard laboratory medium before the onset of the experiments and were never exposed to the media prepared with rotting cactus (see below). We also collected fresh cactus and the juice exudated from naturally occurring rotting cacti of both C. hildmaniannus and P. machrisis. Pieces of fresh cactus were stored at -20°C and the fermenting juices of both maintained in the laboratory by adding 10 g of fresh cactus every 2 weeks.
Since D. gouveai and D. antonietae are difficult to maintain in the laboratory, we could only assess five isofemale lines of each species. All experiments performed were under an LD 12 : 12 h light/dark photoperiod at 25 ± 1°C. Two egg-collecting chambers were set up for each isofemale line. One hundred pairs of sexually mature flies were released into each chamber. Petri dishes containing an egg laying medium were placed in each chamber, removed 12 h later and inspected for the presence of eggs, and incubated for 24 h to allow egg hatching. Batches of 30 first-instar larvae were collected and seeded in culture vials containing 6 mL of a 'seminatural medium' prepared with experimentally rotten pieces of either P. machrisis or C. hildmaniannus.
For the preparation of the cactus media, we followed the methodology described in Fanara et al. (1999). Four replicated vials were initiated for each combination of isofemale line and cactus.
Right wings of four to six males per vial were dissected and mounted on microscope slides. Photographs for each wing were taken using a digital camera mounted on a microscope. Ten landmarks of type I (Bookstein, 1991) at the intersections of veins or where the veins reach the margins of the wing ( Fig. 1) were digitized using TpsDig (Rohlf, 2003a).

MEASURING SIZE AND SHAPE VARIATION
Centroid size was computed as a measure of overall wing size. Centroid size is defined as the square root of the sum of the squared distances of each landmark from the centroid of the configuration (Dryden & Mardia, 1998). Interspecific size differences were tested by analysis of variance (ANOVA) with Species and Cactus as main fixed factors and two random factors: Line (nested in Species) and Replicate (nested in the interaction Line by Cactus).
Sources of intraspecific size variation were analysed by ANOVAs with Cactus (fixed) and Line (random) as main factors and Replicate (random) nested in the Line by Cactus interaction. Besides the mixed nature of the models, the loss of some replicates (summed to the inability of a D. gouveai line to develop in C. hildmaniannus; for a discussion of the biological significance of this observation, see Soto et al., 2007) forced us to use an unbalanced ANOVA design. The model employed, which is specially recommended in cases with missing cells, estimates values for the cells with missing values, using type V sums of squares (Statistica; StatSoft, 2001; see also Hocking, 1996). As a by product of this estimation process, degrees of freedom and the mean squares of the error terms in the ANOVAs are modified.
Shape variation was investigated using the Procrustes technique. This technique is based on the superimposition of all wings and examines differences in the position of the landmarks. Shape coordinates were computed using a least squares Procrustes superimposition method (Bookstein, 1996;Dryden & Mardia, 1998). This method separates shape from size and eliminates variation in scale, position, and orientation (Klingenberg, 2002). The remaining differences in the position of landmarks are due to shape variation (Dryden & Mardia, 1998;Klingenberg & McIntyre, 1998).
The coordinates of superimposed configurations were transformed in shape variables, termed partial warps (Bookstein, 1991;Rohlf, 1996). Partial warps partition shape differences across different geometric scales (Bookstein, 1991) and, because they are nonredundant and geometrically orthogonal, they are amenable for standard multivariate analyses (Rohlf, 1996). The shape space spanned by partial warps and the uniform components of shape variation is termed the tangent space because it lies in a tangent position relative to the curved shape space created by least squares superimposition (Rohlf, 1999). To assess the approximation of tangent space to shape space, the correlation between the Procrustes (shape space) and the Euclidean (tangent space) distances was checked using tpsSmall (Rohlf, 2003b).
Relative warps analysis (Bookstein, 1991;Rohlf, 1993) was used to investigate directions in shape change. This step was performed using a principal component analysis (PCA) of the partial warps and uniform components. As in normal PCA, the relative warps summarize the variation present in a large dimensionality data set in a few major axes. Partial warps and relative warps calculations were performed using tpsRelw (Rohlf, 2003c).
Interspecific shape differences were analysed by means of a multivariate analysis of variance (MANOVA) using the relative warp scores matrix as dependent shape variables (16 relative warps variables) and Species and Cactus as main factors. Intraspecific sources of shape variation were analysed using MANOVAs with Cactus as fixed factor and Line and Replicate (nested in the Line by Cactus interaction) as random factors.

SHAPE VARIATION
The superimposed shape coordinates and the derived warps are not necessarily devoid of allometric change  (Debat et al., 2003). (Debat et al., 2003). To correct the size effect on shape differences (allometric component), we performed inter-and intraspecific multivariate analyses of covariance (MANCOVAs) of the relative warps with centroid size as the continuous covariate. In this way, we tested for non-allometric (size independent) shape differences between the independent factors previously described. We also performed multivariate regressions of the relative warps on centroid size for each species to assess the percentage of shape variation accounted for changes in wing size.

SIZE VARIATION
In the interspecific ANOVA, both main factors, Species and Cactus, were significant ( Table 1). Wings of flies reared in P. machrisis were larger than in C. hildmaniannus (Table 1, Fig. 2) and D. antonietae had larger wings than D. gouveai. The latter showed the smallest centroid size when reared in C. hildmaniannus whereas D. antonietae grown in P. machrisis presented the largest wings (Fig. 2). Random sources of variation were significant in the general ANOVA. Wing size differences among replicates (within lines) may be interpreted as environmental variation. The  significant Line by Cactus interaction points to host dependent genetic variation in wing size within species (Table 1). However, at this point, we cannot determine whether these effects occur in both species or only in one of them.
The intraspecific ANOVAs showed that these two sources of variation (Replicate and Line by Cactus interaction) were significant in both species (Table 1). The Line by Cactus interaction, which may be interpreted as an estimate of the Genotype by Environment interaction (GEI), was significant in both species (Table 1). According to Robertson (1959), a significant GEI may arise due to differences in among-lines variance in separate environments, which is termed a change in scale, and/or deviations of the cross-environment genetic correlation from unity (r GE < 1) (i.e. a change in rank order among lines). The former (change of scale) appears to be responsible for the significant Line by Cactus interaction in D. gouveai because, among lines, variance was greater in C. hildmaniannus than in P. machrisis (Levene's test: F1,182 = 2.14, P = 0.034). In addition, in D. antonietae, variances were significantly different across cacti (Levene's test: F1,207 = 26.83, P < 0.001); however, in this case, changes in the ranking order among isofemale lines across cacti are apparent (Fig. 3). Nevertheless, as observed in the interspecific ANOVA, the significant Cactus effect in both ANOVAs indicated that the rearing hosts had a consistent effect on all lines in both species.

SHAPE VARIATION
The correlation between Procrustes and Euclidean distances was perfect (r = 1.00). Such perfect approximation of shape space by the tangent space made all statistical estimates of shape differences reliable.
The MANOVA detected significant differences in wing shape between species and between flies emerging in different cacti and a significant Species by Cactus interaction (Table 2). Figure 4 plots the mean values of the two main relative warps, which accounted for 46% of total shape variance, for both species reared in both cactus hosts. The shape changes associated with these relative warps (RWs) are depicted as vectors. Both species are clearly differentiated along the RW1 axis, which mainly implied changes in the relative positions of landmarks 2, 3, and 9 (i.e. the relative location of the anterior crossvein and the point where longitudinal vein 5 reaches the wing margin). The intervein region D (IVR D; Birdsall et al., 2000) limited by longitudinal veins 4, 5, and the posterior crossvein (landmarks 4, 5, 8, and 9; Fig. 1) is broader in D. antonietae than in D. gouveai. Moreover, the Species by Cactus interaction detected in the MANOVA is based on the scores matrix for 16 relative warps variables. Even though only the first two relative warps are plotted in Figure 4, the Species by Cactus interaction can still be visualized. In general, the shape of D. antonietae wings emerged in different cacti differed mainly in the scores of RW1, whereas D. gouveai showed remarkable differences among flies reared in different substrates along RW2 axis. Changes across RW1 involves primarily landmarks 2, 3, and 9 shortening or enlarging IVRs C and D. RW2 accounts mainly for variation in the positions of landmarks 6, 4, 5, and 9, although in the latter it is worth to notice that the direction of the vector differs from RW1.
The intraspecific MANOVAs detected significant variation among lines and a significant Line by Cactus interaction in both species (Table 3), suggesting that shape variation is at least partially genetically determined and that phenotypic plasticity for wing shape has also a genetic basis.

SHAPE VARIATION
The multivariate regressions of partial warps on centroid size revealed a significant allometric effect in both species that account for 27.5% and 26.5% of total shape variance in D. antonietae and D. gouveai (P < 0.001 in both cases), respectively. The interspecific MANCOVA showed that shape differences between species, between cactus, and their interaction were also significant after removing size variation (Table 2).
However, the intraspecific non-allometric patterns of shape variation strongly differed between species.  In D. gouveai, all sources of variation remained significant after size correction (Table 3) whereas, in D. antonietae, the effect of the Cactus and the Genotype by Cactus interaction were no longer significant after removing size variation (Table 3). These results suggest that only allometric differences were involved in the morphological differences detected between hosts.

DISCUSSION
Previous studies in the D. buzzatii cluster have addressed the relationship between fitness traits, body size, and cactus hosts in D. buzzatii and D. koepferae. These studies have shown that breeding in the necroses of different cactus hosts may have remarkable consequences on general size, developmental time, survival, and host acceptance (Fanara et al., 1999(Fanara et al., , 2004Fanara & Hasson, 2001;Carreira et al., 2006). Flies reared in the prickly pear Opuntia sulphurea tended to have larger wings, irrespective of the Drosophila species, than those emerged in the columnar Trichocereus terschekii (cardón). However, the most interesting observation is that larval survival in both species was not independent of the rearing substrate. These results suggest that, although both species can live on both types of hosts, there is a certain degree of specialization (Fanara et al., 1999). The authors proposed that cactus chemistry may be involved in host plant specificity as it was also suggested in the guild of Sonoran desert Drosophila (Fogleman, 1982;Fogleman & Abril, 1990). Cactus hosts may differ as breeding sites for the flies in several biologically relevant variables, such as chemical composition (Kircher, 1982;Fogleman & Abril, 1990), diversity in the microflora (yeasts and bacteria) associated with the rotting process (Starmer et al., 1990), rot size, larval density, and rot duration (Hasson et al., 1992). Studies of chemical composition of cactus tissues have shown that columnar cacti are chemically more complex and have lower concentrations of free sugars and, thus, are less nutritious for the flies than Opuntia (Kircher, 1982), hence favouring more specialized, less plastic, flies. Although, in the present study, both cactus hosts are columnars, our results suggest that they may be perceived as different by D. gouveai and D. antonietae. The smaller wing size detected in flies of both species reared in C. hildmaniannus is in line with previous results showing a general poorer performance in this cactus compared with flies reared in P. machrisis (Soto et al., 2007). Estimates of viability for D. gouveai (18%) and D. antonietae (48%) grown in C. hildmaniannus were lower than in P. machrisis (> 80% for both species) (Soto et al., 2007). Considering the available evidence, it may be concluded that C. hildmaniannus may be considered as a relatively stressful environment for the larvae (not only in comparison with P. machrisis) with a detectable impact on adult wing size and general performance. Possible causes of these observations may be related to the presence of toxic compounds in C. hildmaniannus (for similar cases in North American members of the repleta group, see Fogleman, 1982; for the well-  Jones, 1998) and/or to the metabolic activity of bacteria and yeasts living in cactus necroses that may also affect the chemistry of the substrate by increasing or decreasing its toxicity (Fogleman & Danielson, 2001). However, further studies are necessary to identify which features of the cactus chemical composition and/or the microflora associated with the decaying process are involved in the contrasting responses of flies to the cactus hosts offered as feeding resources in the present study.
Wing morphology also showed a plastic response in both species, particularly for wing size, but was less evident for shape. Although the number of lines analysed was not large, the Line by Cactus interaction, which may be construed as an estimate of the genotype by environment interaction, detected in both species suggests the presence of genetic variance for wing plasticity. The relationship between size and shape variation of the wing and the possibility of independence between the two components of wing morphology has been the subject of much discussion (Weber, 1990(Weber, , 1992Bitner-Mathé & Klaczko, 1999a, b;Weber et al., 2001;Hoffmann & Shirrifs, 2002;Carreira et al., 2006). Shape changes accompanying size increments (allometric changes) are well documented in insect wings (Stern & Emlen, 1999). In our case, wing shape differences between flies reared in both cactus hosts (i.e. phenotypic plasticity) disappeared after removing the effect of size in D. antonietae, indicating that wing shape variation in this species can be mainly accounted for by the allometric relationship between size and shape. On the other hand, phenotypic plasticity in wing shape showed two components in D. gouveai: one that can be explained by the allometric relationship between size and shape and a non-allometric component. This is an interesting observation that agrees with what would be expected. When moving a species such as D. antonietae from its own environment (although stressful) to a better environment, the effect is preferably detectable on size because developmental noises are not expected. By contrast, when D. gouveai is moved from its own environment to the stressful one, we would expect effects on both size and shape (non-allometric). All in all, interspecific differences in the patterning of the plastic response have probably evolved because these species shared their last common ancestor 3.3 Mya (estimated from DNA sequence data reported in Manfrin et al., 2001). Moreover, this is the first study to report host-related plasticity in wing shape in the D. buzzatii cluster because the non-allometric component of wing shape did not contribute to variation in wing morphology in other species of the D. buzzatii cluster (Carreira et al., 2006). A significant proportion of variation in wing size and shape can be accounted for by differences among lines and the Line by Cactus interaction in both D. gouveai and D. antonietae. Similar results were obtained for viability and developmental time in the same set of isofemale lines (Soto et al., 2007). These results indicate a potential role of environmental heterogeneity (as represented by different host plants) in the maintenance of genetic variation in wing morphology. This function of environmental heterogeneity appears to be widespread, at least in the D. buzzatii cluster (Carreira et al., 2006;Fanara et al., 2006), indicating that host plant use played a significant role in the evolutionary history of the D. buzzatii cluster.
Previous studies have shown that wing morphology provides useful traits for species discrimination in the D. buzzatii cluster (Moraes et al., 2004a). Moreover, it was concluded that wing morphology is a phylogenetically informative trait because phenetic relationships among species, based on wing traits, are compatible with phylogenetic trees inferred using molecular data. Our study shows that divergence in wing morphology between D. gouveai and D. antonietae is greatest in IVR D, which includes the posterior crossvein (Fig. 1, landmarks 4, 5, 8, and 9). Interestingly, this portion of the wing is also a species-specific trait, useful for the identification of closely-related species such as D. mercatorum and D. paranaensis (Moraes et al., 2004b) and D. buzzatii and D koepferae (I. M. Soto, unpubl. data). Moreover, IVRD explains the major part of the plastic response observed in D. antonietae and D. gouveai in response to the cactus rearing media.
The present study may also contribute to unveil certain aspects of wing development and evolution that appear to be elusive. The primordium of the Drosophila wing is divided during early development into distinct portions (García-Bellido & De Celis, 1992;Resino, Salama-Cohen & Garcia-Bellido, 2002) and studies of wing morphology have shown that the major compartments of the wing behave as independent units of selection contributing differentially to adaptation (Cavicchi et al., 1985). These observations have led to the proposal that intervein regions may represent morphological modules that are under different genetic control, implying that growth properties in different portions of the wing may be spatially coordinated Pezzoli et al., 1997). In the present study, suggestively, the same region of the wing that appears to be more prone to interspecific divergence was also involved in the major shape changes involved in the plastic response experienced by both species exposed to different host cacti (Fig. 4). This observation is in line with previous studies in D. melanogaster that reported a more plastic response of posterior wing compartments to environmental variation of different source (Cavicchi et al., 1985;Guerra et al., 1997;Pezzoli et al., 1997). Although it may be premature to advance an explanation for these observations, our results suggest, at the very least, that the evolution of different parts of the wing is under the influence of different constraining factors.