Butterflies possess attributes that are sensitive to gradual environmental changes. Recently, the effects of environmental factors on the shapes of organisms, as well as the interactions of these elements, have been extensively examined, i.e., effects of seasonal changes on the colors of butterfly wings, and effects of landscape structure on butterfly distribution and morphology. However, few studies have dealt with variations in butterfly shapes in response to varying environmental conditions. Here we aimed to determine how body size and shape variations in butterflies are correlated to environmental heterogeneity. We used geometric morphometrics to quantify Pieris rapae wing shape variations. Results showed that forewing and hind wing sizes were significantly different among the 15 populations. P. rapae individuals with larger wing sizes were mainly distributed in mountainous areas, whereas those with smaller-sized wings were found on the plains. Canonical variate analysis was employed to examine the patterns of variation in wing shapes among and within the populations. Significant differences in shape were revealed in the forewings and the hind wings of P. rapae populations. All populations were divided into 2 groups on the first canonical variate axis (CV1), which followed the Qinling Mountains as an important boundary between the Palearctic and Oriental Realms in zoogeographical division of the world. The unweighted pair group method with arithmetic mean (UPGMA) clustered the 15 populations into 4 groups by forewing and hind wing shape in response to the 4 environment types in Qinling Mountains. We suggest that wing shapes of P.rapae are sensitive to environmental heterogeneity. The isolating effect of the Qinling Mountains on P. rapae population interactions was apparent.
Understanding how environmental heterogeneity affects the phenotypic patterns of organisms is a major focus of evolutionary ecology (Monaghan 2008; Fischer et al. 2010; Moczek 2010; Vargas et al. 2010). Under certain environmental conditions, changes in an organism's phenotype can increase its fitness, and thus organisms exhibit the capacity to adjust their phenotype to match prevailing local conditions (Merckx & Van Dyck 2006; Monaghan 2008; Otaki et al. 2010). Recently, the effects of environmental factors on the shapes of organisms, as well as the interactions of these elements, have been extensively examined (Beldade & Brakefield 2002; Prieto & Dahners 2009). These effects on organisms include the following: food resource effects on the horns of dung beetles (Pfennig et al. 2010), effects of seasonal changes on colors of butterfly wings (Daniels et al. 2012), photoperiod and temperature effects on the body size of grasshoppers (Harris et al. 2012), as well as the effect of landscape structure on butterfly distribution (Vandewoestijne & Van Dyck 2011) and on butterfly morphology (DeVries et al. 2010). However, few studies have linked animal shape variations to changing environmental conditions. The reasons for this lack of studies include the complexity and diversity of morphological variations. Geometric morphometrics has developed and matured sufficiently to support this important branch of morphological research.
Geometric morphometrics has been applied in the quantitative analyses of shape variation. Shape variables are computed and regressed onto geographical coordinates and environmental variables by both linear and curvilinear models (Cardini et al. 2007). In the present study, we use geometric morphometrics to examine the variation in size and shape of butterfly wings in diverse environments. The purpose is to understand the relative roles of shape variation and environmental effects, as well as the interactions between them, in shaping geographical population patterns.
Butterflies, including our model species, Pieris rapae L. (Lepidoptera: Pieridae), have been thoroughly examined for their developmental and phenotypic variations in life-history traits and adult morphology across space and time(Fric et al. 2006; Breuker et al. 2007; Gibbs et al. 2011). Butterflies are also known to be highly sensitive to climate (Dennis 1993). In large geographical areas, vegetation and climate often vary greatly, especially those regions blocked by large mountain ranges. Changes in temperature affect all aspects of butterfly life history, including their distribution and abundance. Changes in rainfall levels can indirectly influence butterfly larvae through changes in host plant quality (Roy et al. 2001; Morecroft et al. 2002). Butterflies may alternatively allow for plasticity in their phenotype, linked to environmental variation. P. rapae is a sexually dimorphic species, and recent research has demonstrated that the female of this species is more likely to exhibit much morphological variability (Stoehr & Goux 2008; Snell-Rood & Papaj 2009). Wing patterns in butterflies are not only visually stunning examples of the evolutionary process, they are also emerging as exceptional model systems linking developmental and genetic processes, generating morphological variations with ecological and evolutionary processes, and thereby molding variations in natural populations(McMillan et al. 2002; Brakefleld 2006; Brakefleld et al. 2007). Environmental stress can affect the genome and expression of genetic variation at the butterfly phenotypic level, especially the wings. Wing shape exhibits high degrees of variation at inter- and intraspeciflc levels (Kingsolver 2000, 2004; Talloen et al. 2009, Dincă et al. 2011). Minor variation in butterflies can be analyzed quantitatively at the species and subspecies levels (Prieto et al. 2009; Miguel et al. 2011; Kitthawee & Rungsri 2011).
In this study, we selected the wing shape and venation of P. rapae for analysis. These 2-dimensional structures are stable, and the veins are clearly visible (Betts & Wootton 1988; Nygren et al. 2008). Homologous landmarks on the nodes of wing venation and wing edges help describe the differences in wing shape among different geographical populations. Here, we considered the hypothesis that environmental forces influence the sizes and shapes of the wings of P. rapae as our main point for evaluating shape variations of butterflies in diverse environments.
Materials and Methods
STUDY AREA AND SAMPLING
All P. rapae samples were collected from the Qinling Mountains and adjacent regions. The Qinling Mountains form a natural boundary between the north and the south of the country. The northern side of the range is prone to temperate weather, whereas the southern side has a subtropical climate with a rich, fertile landscape supporting abundant wildlife and vegetation. The Qinling Mountains are part of the boundary between the Palearctic Realm and the Oriental Realm, and these tall mountains block the interflow of species. Hence, the Qinling Mountains are an ideal place to study diverse environments and shape variations in organisms.
We collected specimens in the Qinling Mountains and adjacent regions from Jun to Aug in 2008. A total of 15 geographical populations (Fig. 1) were identified from north of Shaanxi, the Guanzhong Plain, south of Shaanxi, and around the branch of the Qinling Mountains in the Henan areas. We used a geographic information system (GIS) and its distance measuring tools to calculate the geographic distance between each of populations. Among sites of these populations, the longest straight line distance was 550 km, whereas the shortest was 25 km. The following conventions were used to classify the P. rapae populations with respect to types of environment (Table 1): 1 — population from hilly and artificial vegetation environments with a warm temperate semi-arid climate; 2 — population from the transition zone between a warm temperate semi-arid climate and a warm temperate sub-humid climate; 3, 4, 5, and 6 — populations from areas with a warm temperate sub-humid climate; 8, 9, 10, 11 — populations from mountainous, and natural vegetation environments with a cool subtropical humid climate; 7 — population from the transition zone between a warm temperate sub-humid climate and a cool subtropical humid climate; 12, 13, 14, 15 — populations from the plains in a transition zone between a warm temperate sub-humid climate and a warm temperate humid climate.
Fig. 1.
Distribution map of the P. rapae populations studied and the integrated physical regionalization (diverse environments) in the Qinling Mountains and adjacent regions. Note: The numbers represent the IDs of the populations; the circles and groups represent the populations divided by the cluster analysis from Fig. 6.
DATA ACQUISITION
Wing size and shape variations were examined and recorded from at least 20 female individuals per location by the landmark based geometric morphometric method (Bookstein 1991; Rohlf et al. 1996; Adams et al. 2004) for a total of 300 specimens of P. rapae from the 15 geographical populations. Images of the right forewing and right hind wing of each female specimen were captured using a Sony DSCH5 camera attached to the copy stand, with a fixed focus and the same camera angle and magnification ratio for all specimen images captured. A total of 14 landmarks on the forewing and 12 landmarks on the hind wing positioned at vein intersections or terminations (Fig. 2) were collected and digitized using TpsDig 2.10 (Rohlf 2006). These landmarks were used to correspond to x, y coordinates in a Cartesian space (Adams et al. 2004).
MORPHOMETRIC AND STATISTICAL ANALYSES
The specimen's wing size (measurement unit: mm) was calculated based on the centroid size (CS; the square root of the sum of squared distances between each landmark and the wing centroid), an isometric estimator of size (Zelditch et al. 2004). The differences in centroid size (CS) among populations were analyzed by one-way analysis of variance (ANOVA with post hoc Tukey's HSD test). To examine wing shape variation, digitized landmark data were subjected to generalized procrustean superimpositions to standardize the size of the landmark configurations and eliminate differences due to translation and rotation (Rohlf & Slice 1990). Major shape changes in projected lateral view were illustrated using thin-plate spline analysis (Bookstein 1989). The resulting weight matrix (Rohlf 1993) was then used to explore shape change by means of a multivariate canonical variate analysis (CVA). Furthermore, the visual representation of shape differences described by canonical variates was produced by regressing the shapes (the weighted matrix of the partial warp scores) onto the specimen scores on the first 2 canonical vectors (Zelditch et al. 2004; Klingenberg et al. 2013). This permitted the splines of the shape change to be associated with positive and negative values of a canonical vector.
Table 1.
Sampling localities, population identifiers (ID, used in Figs. 1, 3, 4, 5 and 6 and in the main text) employed in this study.
The morphometric analyses were conducted using the IMP software package for geometric morphometrics (Rohlf 2006), and CVA was performed in PAST ver. 1.75 (Hammer & Harper 2007). The relationship between the morphometric shape characteristics of the geographical populations was illustrated by cluster analysis using the unweighted pair group method with an arithmetic mean (UPGMA) based on Euclidean distances between mean shapes (computed from partial warp scores between pairwise population consensus configurations). Cluster analysis was performed using standard algorithms of NTSYS-pc v.2.10e program (Rohlf 2002). The above- mentioned consensus configuration of each population was performed using the IMP software (Sheets 2000).
Results
WING SIZE DIFFERENCES WITHIN 15 PIERIS RAPAE POPULATIONS
The Shapiro-Wilk test revealed a normal distribution in all populations in size of the forewing and hind wing (P > 0.05). One-way ANOVA of the mean forewing CS (Fig. 3) showed significant differences between inter-population variations (F(14, 298) = 5.98, P = 0.001). Forewing size differed significantly among 3 groups of populations: the character CS varied significantly between populations 3 and 4 (larger wing size), as well as among populations 12, 13, 14, and 15 (smaller wing size) (P < 0.001); the character CS varied significantly between populations 8 and 9 (larger size) and among populations 12, 13, 14, and 15 (smaller size) (P < 0.001); the character CS varied significantly among populations 3 and 4 (larger wing size) and population 5 (smaller wing size) (P < 0.001). One-way ANOVA of the mean hind wing CS (Fig. 3) showed significant differences between inter-population variations (F14, 298 = 7.07, P = 0.001). The size variance of the hind wing and the forewing showed a high level of similarity by the same analysis method, and this observation indicated that the character CS of populations 3, 4, 8, and 9 (larger wing size) was significantly different from that of populations 12, 13, 14, and 15 (smaller wing size) (P < 0.001). The character CS varied significantly among populations 3 and 4 (larger wing size) and population 5 (smaller wing size) (P < 0.001).
WING SHAPE VARIATIONS OF THE 15 PIERIS RAPAE POPULATIONS
The results of the geometric morphometric analysis of the wing shape were visualized by CVA and the thin-plate spline analysis (Figs. 4 and 5). The CVA of shape variability among the 15 populations clearly showed that the differences in the forewing were highly significant among populations (Wilks' Λ = 0.04, F392, 3289 = 2.41, P < 0.001); however, the populations overlapped in all scatter plots. The first 2 axes of the CVA exhibited 46.68% and 14.88% of the forewing variation. Meanwhile, the hind wing shape variability among populations indicated a signiflcant difference (Wilks' Λ = 0.07, F336, 3224 = 2.27, P <0.001). The first 2 axes of the CVA exhibited 50.41% and 15.49% of the hind wing variation. We selected the average centroid distributions of each population to explain the shape variance among the populations because of the large overlap in the populations. In general, the shape of the forewing and the hind wing of the populations on the first axis (CV1) could be used to divide the 15 populations into 2 groups, which is consistent with the characteristics of the Qinling Mountains as the boundary between northern and southern China. Populations 8, 9, 10, 11,12, 13, 14 and 15 were mainly distributed on the CV1 positive axis as the South group, whereas populations 1, 2, 3, 4, 5 , 6 and 7 were mainly distributed on the CV1 negative axis as the North group. Thin-plate spline analysis showed that forewing shape deformation was mainly derived from the discoidal cell and the R vein (landmarks 2, 3, 4 and 5) and between the M2 and M3 veins (landmarks 10 and 11). Hind wing shape deformation was mainly derived from the discoidal cells (landmarks 1, 2, 3, 4 and 5) and between the M2 vein and M3 vein (landmarks 8 and 9) (Fig. 2, Figs. 4 and 5).
The forewing and hind wing shapes of the populations on the second axis (CV2) were not clearly distinguished (Figs. 4 and 5). Thin-plate spline analysis showed that forewing shape deformation was mainly derived from the base of the wing (landmarks 1 and 6) and between the R3 and M1 veins. Hind wing shape deformation was mainly derived from the base of the wing (landmarks 1 and 6) and between the Cu1 and Cu2 veins (Fig. 2, Figs. 4 and 5).
The aforementioned results show that forewing shape variance was highly similar to that of the hind wing (Figs. 4 and 5). The 15 populations could be divided into 2 groups based on the CV1 axis. The populations south of the Qinling Mountains were distributed on the positive axis whereas the populations north of the mountain were distributed on the negative axis. Thus the first canonical variate axis corresponds to the Qinling Mountains as an important boundary between the Palearctic Realm and the Oriental Realm in the zoogeographical division of the world. The populations on the CV2 axis were not clearly distinguished, but the shape variance shows that the populations in the same environments shared a similar wing shape character. Examples include populations 9 and 10, populations 2 and 3, as well as populations 10 and 11.
WING SHAPE RELATIONSHIPS AMONG THE 15 PIERIS RAPAE POPULATIONS
UPGMA cluster analyses were used to evaluate the shape relationships of the 15 populations, whereas cluster data from Euclidean distances were computed for the partial warp scores between mean shapes of each 15 population (consensus configurations). The results (Fig. 6) revealed that the forewings of the 15 populations could be clustered into 4 groups with a linkage distance of 0.0027: Group A, population 1 as one branch; Group B, populations 10, 11, 12, 13, 14 and 15 as one branch; Group C, populations 2, 3, 4, 5, and 6 as one branch; Group D, populations 8 and 9 as one branch; population 7 isolated as one branch, but very closely related to group C. The hind wings of the 15 populations could be clustered to 4 groups with a linkage distance of 0.0038. The cluster analysis results were similar to those of the forewing except that population 7 was clustered with Group C as one branch. The forewing and hind wing cluster analysis results were consistent with 4 kinds of environmental types in the Qinling Mountains and adjacent regions (Fig. 1). Group A (population 1) was located in the warm temperate semi-arid climate as well as in hilly and artificial vegetation environments; Group B (populations 10, 11, 12, 13, 14 and 15) was located in the transition zone between warm temperate subhumid and humid climates of the plains; Group C (populations 2, 3, 4, 5 and 6) was located in the warm temperate sub-humid climate of the plains; Group D (populations 8 and 9) was located in cool subtropical humid climate, mountainous region with natural vegetation environments.
VARIATIONS IN WING SHAPE AND GEOGRAPHICAL DISTANCE
We selected the populations 2, 3, 4, 5, 6, 13, 14 and 15 to analyze the relationship between wing shape variance and geographical distance. The locations of all aforementioned populations had the same climatic conditions and were not isolated by mountain ranges (Fig. 1). The correlation between the geographical distance between populations and the Euclidean distance between wing shapes for each forewing and hind wing was analyzed. The geographical distance in each population was an independent variable. The Euclidean distance (shape differentiation) in each of population was a dependent variable (Fig. 7). The results showed that the Euclidean distances for the forewing and the hind wing had respectively significant positive correlations with the geographical distances (P < 0.001, r = 0.58, forewing; P < 0.001, r = 0.78, hind wing). This means of the shape variances of P. rapae populations increased with increases in geographical distance.
Discussion
The butterfly P. rapae exhibits changes in wing size (CS) across environments in the Qinling Mountains (Fig. 3, Table 1). The larger-winged populations were shown to be mainly distributed in mountainous areas, whereas those with smaller-wings were mainly distributed on the plains and suburb of cities. This wing size difference is probably influ enced by environmental factors such as food, terrain, and city effects. Vargas (2010) suggested that the wing sizes of insects are associated with the growth of the larva. Larval growth is directly affected by temperature, and to a greater degree, by humidity. Insects can develop larger wings in environments with increased humidity. Moreover, insects often encounter barriers in mountainous areas that are absent on the plains. The large wings may be favorable for finding mates, food sources and adapting to specific environments (Prieto & Dahners 2009). The smaller P. rapae wings evident in the suburbs may be due to city effects. Kingsolver & Huey (2008) suggested that P. rapae often exhibited smaller sizes in high-temperature environments. The city effects lead to increased temperatures and reduced humidity, which may be the reason for the smaller wing size of the P. rapae population near the city. This conclusion agrees with the previous research of Schoville (2013). However, the butterflies' sizes differ somewhat between plants and between seasons. Host effects are seen as an important factor influencing the size of larvae, which may result in size differences in adults (Dennis et al. 2005; Friberg & Wiklund 2009). Whether the wing size of P. rapae clinal variation is caused by city effects or host effects, the situation seems complicated and further research will be necessary to determine causation.
Fig. 3.
Boxplot of P. rapae centroid size (CS) with the mean, standard error, and standard deviation illustrating variations in wing size across geographical populations.
Fig. 4.
Forewing shape variation (CV1: 46.68%; CV2: 14.88%). The colored circles in the image above represent the average discrete point centers of populations; the number is the population ID. Thin-plate spline analysis results are shown by colored grid, which represents wing shape deformation. The numbers on the grid are landmarks of wings. Blue colored notes contraction between landmarks, and red color indicates expansion between landmarks. The North group and South group correspond to the boundary of Qinling Mountains as the boundary between northern and southern China.
Fig. 5.
Hind wing shape variation (CV1: 50.41%; CV2: 15.49%). The colored circles in the image above represent the average discrete point center of populations; the number is the population ID. Thin-plate spline analysis results are shown by colored grid, which represents wing shape deformation. The numbers on the grid are landmarks of wings. Blue color denotes contraction between landmarks, and red color indicates expansion between landmarks. The North group and South group correspond to the boundary of Qinling Mountains as the boundary between northern and southern China.
Environmental diversity in Qinling Mountains contributes to P. rapae wing shape variation. The dramatic environmental variation from north of the Qinling mountains to the south of these mountains allowed an in-depth exploration of the wing shape variation in P. rapae populations. The CVA of wing shape variability among the 15 P. rapae populations can be divided into a north group and a south group by phenotypic differences plotted on the CV1 axis (Fig. 4, forewing and Fig. 5, hind wing). The groups are consistent with the characteristics of the Qinling Mountains being the boundary between northern and southern China: the northern group belongs to the Palearctic Realm and the southern group belongs to the Oriental Realm.
This demonstrates the role of the Qinling Mountains as a barrier that has driven intraspecific fragmentation. It was highly unexpected that wing-shape variation in populations of P. Rapae would be so perfectly in line with the existing geographical barriers. Thin-plate spline analysis showed that both the forewing and the hind wing of the northern group have a larger discoidal cell and smaller medius area, and the southern group has a smaller discoidal cell and larger medius area. The deformation of the P. rapae forewing and hind wing mainly occurs at the intersection between the discoidal cell and the medius area, possibly affecting P. rapae flight ability. For the different life histories in diverse environments, selection should act on wing shape to optimize flight capabilities (Sunada 1993).Some studies show the insect wings are stiffer at the base than at the margins. Greater flexural strength is observed in the discoidal cell and medius area (Combes & Daniel 2003a, 2003b). Here we predict that the discoidal cell and the medius area may be the mediators of force during P. rapae flight.
Cluster analysis (Fig. 6) divided the wing shapes of the 15 populations of P. rapae into 4 groups. The groups are consistent with the main environments types in the Qinling Mountains (Fig. 1). Such results suggest there is a close relationship between shape variation and diverse environments within the range of distribution of P. rapae. Wing shape variability is an important means by which butterflies cope with environmental variations (Brakefield & Reitsma 1991; Schlichting & Pigliucci 1998).Microhabitat selection, such as larval development in some species of butterflies has been shown to be temperature- and humidity-sensitive (Merckx et al. 2003). This study showed that the various P. rapae populations distributed around the Qinling Mountains encounter a complex geography and landforms, special geographical conditions and various climates. Individuals in P. Rapae populations exhibited diversified morphological shapes in these different conditions. Previous studies of the giant panda Ailuropoda melanoleuca (Zhang et al. 2002), birds (Lei et al. 2003) and amphibians (Chen et al. 2008) in the Qinling Mountains revealed that environmental factors also contribute to differentiation of shapes and genetic structures of these animals. This demonstrates that the isolation effects of the Qinling Mountains on population interactions are universal and enable organisms to adapt to their niches and evolve independently.
Fig. 6.
UPGMA tree of P. rapae forewing and hind wing among different populations, based on Euclidian distances between mean wing shapes. The cluster numbers are population IDs (see Table 1).The groups are divided by Euclidian distances, i.e., the forewing divided by a linkage distance at 0.0027 and the hind wing by a linkage distance at 0.0038.
In this study, the wing shape variation of P. rapae was highly correlated with geographical distance: the smaller geographical distances among populations with smaller shape differences (Fig. 7).These results suggest that the geographically proximate populations of P. rapae exhibit a little shape differentiation, but as the distance among the populations increases, shape differentiation among populations also increases. In general, populations that are farther apart have a smaller chance of gene exchange. Thus, shape change usually occurs in insect populations because of isolation by distance (Kingsolver et al. 2007). Pieris rapae are sensitive to habitat fragmentation (Schoville et al. 2013).Their distribution and interactions among their various populations may be related to population-specific habitat use and requirements and thus the scale at which they perceive the structure of the landscape.
Acknowledgments
We are much obliged to the 2 anonymous reviewers for their constructive suggestions on improving this manuscript. We sincerely thank John Richard Schrock of Emporia State University, USA, for help in improving the English of this manuscript.
This work was supported in part by the National Natural Science Foundation of China under Grant No. 31372250 and No.31402006, was also supported by the Fundamental Research Funds for the Central Universities under Grant No. GK201302053.
