Patriline Differences Reveal Genetic Influence on Forewing Size and Shape in a Yellowjacket Wasp (Hymenoptera: Vespidae: Vespula flavopilosa Jacobson, 1978)

The wing venation is frequently used as a morphological marker to distinguish biological groups among insects. With geometric morphometrics, minute shape differences can be detected between closely related species or populations, making this technique useful for taxonomy. However, the direct influence of genetic differences on wing morphology has not been explored within colonies of social insects. Here, we show that the father’s genotype has a direct effect on wing morphology in colonies of social wasps. Using geometric morphometrics on the venation pattern, we found significant differences in wing size and shape between patrilines of yellowjackets, taking allometry and measurement error into account. The genetic influence on wing size accounted for a small part of the overall size variation, but venation shape was highly structured by the differences between patrilines. Overall, our results showed a strong genetic influence on wing morphology likely acting at multiple levels of venation pattern development. This confirmed the pertinence of this marker for taxonomic purposes and suggests this phenotype as a potentially useful marker for phylogenies. This also raises doubts about the strength of selective pressures on this phenotype, which highlights the need to understand better the role of wing venation shape in insect flight.


Introduction
Morphology is the result of complex interactions between genetic and environmental influences on the organism during development [1]. Being the most obvious interface between the organism's genetic background and the environment, morphological characters have long been studied to classify organisms and to understand their evolution. Among insects, wing venation is one of the most studied structures due to its strong taxonomic signal [2,3], its frequent occurrence in the fossil record [4,5], and its functional importance [6,7]. With the development of geometric morphometric methods, more and more studies use wing venation shape as a morphological marker to distinguish groups, including species of bees and wasps [3,5,8]. Wing venation shapes differ between closely-related species, but also populations, sex and castes within populations [9][10][11][12]. Monteiro and collaborators [13] even found significant variation of this marker between colonies of honeybees, suggesting the possibility of intraspecific genetic effects on venation. However, while they inferred an estimate of the wing shape heritability, their data enabled no direct observation of the influence of the genotype on this phenotype.
The wing is known to be influenced by external factors such as temperature, diet and cell size [12,14,15]. Furthermore, wing shape variation tends to decrease as genetic distance decreases across taxa [3]. The relative proportion of the shape variation caused by these external factors is thus expected to increase when genetic variation lowers. The observed differences in wing shapes could thus be caused by environmental variation as much as genetic differences, which could be a problem when taxonomic groups are determined with this marker [16,17]. Genetic effects on insect wing shapes have been detected in Drosophila through mutation accumulation [18]. However, these results were restricted to organisms reared in artificial conditions, and recent studies have shown that the genetic process of venation formation is different in Hymenoptera [19]. In order to assess the importance of genetic variation on wing shape in the Hymenoptera, we tested for the first time whether patrilines could be distinguished on the basis of their wing venation.
The social Hymenoptera are particularly useful subjects for the study of the genetic basis of morphology [20][21][22]. In species like honeybees and yellowjackets, nests contain numerous adults derived from a single mother, the queen. For species with queens that mate with multiple males, identifying full-sibling groups of daughters (patrilines) is straightforward, once the queen's genotype is known. Furthermore, the father's genetic influence is identical over its entire patriline because male Hymenoptera are haploids. The existence of multiple genetic lines within the common environment of the yellowjacket nest allows the measurement of the relative roles of genotype and environment in explaining phenotypic variation.
In this paper, we explore the genetic effects on wing size and shape of the yellowjacket wasps Vespula flavopilosa (Jacobson, 1978). We first determined paternity for workers from two colonies using microsatellite markers before testing whether wing size and shape differed between these patrilines.

Sampling
The two mature colonies of V. flavopilosa were collected in Ithaca (NY, USA, 42°27'10.25"N 76°28'36.78"W) in late summer of 2013 (Colony 1: August 23 rd ; Colony 2: September 10 th ). The nests were collected on private land with verbal authorization from the land owner. No specific permits were required for the collection of this abundant, non-endangered species. Colonies were anesthetized overnight using CO2 prior to the collection. Adult workers were then stored at -20°C until DNA was extracted. These colonies were chosen from among 10 colonies, studied in a previous analysis of the species [23], based on the number of available workers and the number of observed patrilines. by employing simple dye-labeled primers (Applied Biosystems) and a 3-primer method [27], as well as fragment size for those labeled with the same dye. Each 10 μl PCR reaction included 1 μl of extracted DNA, 5 μl Qiagen master mix (Qiagen Type-It Microsatellite Kit, Qiagen Inc.), 0.2 μl of each reverse primer, 0.2 μl (dye-labeled) or 0.1 μl (3-primer labeled) of each forward primer, 0.15 μl FAM-labeled 3-primer tag for each 3-primer-labeled primer pair, and water. PCR reaction conditions were 95°C for 15 minutes, 35 cycles of 95°C for 30 seconds, 50°C for 90 seconds, 72°C for 60 seconds, followed by 60°C for 30 minutes. The fragment analysis was performed on an ABI-3730xl sequencer using 0.5 μl PCR product combined with 15 μl HiDi Formamide and 0.15 μl LIZ 500 internal size standard (Applied Biosystems). Allele sizes were called using GeneMarker (SoftGenetics LLC) and checked twice by eye.

Assigning paternity
We used Colony2 v2.0.4.1 [28] to find the maximum likelihood configuration of paternity assignments for all genotyped workers. Workers that failed to amplify at more than one locus were excluded from the analysis. We included 116 genotyped workers from Colony 1 and 112 genotyped workers from Colony 2, as well as 157 workers from the eight additional colonies reported in the previous study [23]. In the few cases where Colony2 initially assigned a worker to a matriline from a different colony, the anomalous worker genotype was checked against the genotype of the queen from that colony. If she shared an allele with the queen at all loci, a maternal sibship constraint was entered into Colony2 containing all workers in that colony that did not differ from the queen for both alleles at any locus. In all cases, the subsequent run of the likelihood analysis assigned the worker to an additional patriline from that colony.
The probability that two randomly selected males shared the same multi-locus genotype (non-detection error) was assessed using the formula where q ij is the frequency of the ith allele at the jth locus [29]. Allele frequencies were calculated and adjusted by Colony2.

Morphometric measurements
The right forewing of every specimen was dissected and mounted in a water droplet between microscopic slides. In order to obtain repeatable pictures of wings that were as coplanar as possible, we removed the thick part at the base of the wings beforehand. Pictures were taken with a Visionary Digital imaging system, using Infinity lenses and a Canon 60D camera with a constant magnification. Forewing size and shapes were assessed using 19 2D landmarks recorded with TPSDig2 software [30] (Fig 1). Landmarks of all specimens were superimposed once, using a generalized Procrustes analysis (GPA; [31]). This analysis extracts information of location, orientation and size from the raw landmark coordinates measured on the pictures in order to retain only the geometric shapes. Aligned landmark coordinates were then projected into the linear tangent space [32]. Shape variables were assessed as the scores of the 34 Principal Components (PC) with eigenvalues greater than zero from a Principal Component Analysis (PCA) on these tangent coordinates. In this multidimensional space, each point represented a potential wing shape and each vector depicted a shape change.
The forewing size was estimated by the log-transformed centroid size from the Procrustes analysis [33].

Measurement error
As the observed shape variation was very low, we verified that the variation due to measurement error was lower than the biological variation. Measurement error was divided between two potential sources of variation: the preparation of the wings and the accuracy of the digitization. In order to compare the scale of this variation relative to the biological variation, the wing of 39 specimens from Colony 1 were prepared twice and landmarks were digitized twice for each picture. The scale of the measurement errors was assessed by comparing the mean squares of wing shapes variation [34]. Mean squares were computed from a Procrustes Analysis of Variance (Procrustes ANOVA) using hierarchical system of four nested factors: patriline, individual, wing mounting and landmark digitization.

Analyses
First, the variation in wing shape between colonies and between patrilines was observed on the whole sample. The wing size difference between colonies and between patrilines was tested using an analysis of variance (ANOVA) with two nested factors: colony and patriline. Allometry, the effect of size on shape, and wing shape differences between colonies and between patrilines within colonies were tested using a multiple analysis of covariance (MANCOVA). The actual difference in wing shape between patrilines was further estimated by the rate of correct reassignments using canonical variates with leave-one-out cross-validation. Reassignments were estimated using 2 to 34 PCs because canonical variates may be influenced by the dimensionality [35].
Because several patrilines had low sample sizes, further analyses were restricted to balanced subsamples (34 individuals) drawn from the two main patrilines in each colony. Wing size differences between these patrilines were tested in each colony using Student's t-tests. The influence of patrilines on wing shape was tested on both shape data and allometry-corrected data using the two-sample Hotelling's T² test. Allometry-corrected data were assessed as the residuals of colony shape data from the common allometric component. The common allometric component was estimated for each colony by a multivariate regression of pooled patriline shape variables on size. Shape differences between colonies and between the main pair of patrilines of each colony were also compared in terms of angles between the vectors within the tangent space. If the differences between these three groups were similar, the three vectors describing these shape changes should have a similar direction. We thus tested whether the angles between these vectors of group differences were significantly lower than angles between pairs of shapes randomly drawn from the sample.

Results
We detected four patrilines in each colony. Non-detection error was low (1.5%), suggesting a low probability of mistakenly identifying workers from two distinct patrilines as the offspring of a single male. Patriline frequencies varied from 14.66% to 37.07% in Colony 1 and from 5.36% to 42.86% in Colony 2. Out of the 228 successfully genotyped specimens, eight had the wing too damaged for proper measurements. A total of 220 wings were thus measured for size and shape (Table 1).

Shape variation
Measurement error was more than a hundred times lower than the differences in wing shape between patrilines ( Table 2). The variation between individuals within a patriline was five times higher than the error induced by the wing preparation, and more than 292 times higher than the error related to landmark digitization. These results confirmed that the observed biological variation of the wing shape could not be explained by measurement error.
The first two principal components (PCs) of the shape variation of the whole sample separated the two colonies and revealed a partial structure of the different patrilines within colonies (Fig 2). These PCs accounted for 33.5% and 13.3% of the total shape variation respectively. Additional differences between pairs of patrilines were observed on several of the other PCs. The MANCOVA showed that allometry-the influence of size on the wing shape-was significant, as were the shape differences between colonies and among patrilines ( Table 3). The allometry was significantly different between colonies, but not among patrilines within colonies. The average procrustes distance between patrilines' mean shapes was only slightly higher than the average procrustes distance between an individual shape and the mean shape of its patriline (0.01422 ± 0.00399 between patriline mean shapes; 0.00951 ± 0.00212 within patrilines). However, the canonical variates confirmed that different patrilines had different wing shapes: the rate of individuals that could be attributed to their correct patriline on the basis of their wing shape was high. When using cross-validation, it ranged from 65% when using 2 PCs-46.6% of the total shape variation-to 99.09% with 19 PCs-96.53% of the total variation (S1 File). The wing shapes of different patrilines thus occupied overlapping, but largely distinct, regions of the shape space. Within each colony, the two main patrilines were also readily discriminated on the first PCs of colony-based PCAs. The Hotelling T² tests confirmed that the wing shapes of individuals Two first Principal Components (PC) of the shape variation of the wing venation patterns. The region occupied by each patriline was indicated in grey. The shape changes described by each PC were amplified threefold to permit visualisation. Percentages of the axes indicate the proportion of the total shape variation explained by each PC. The apparent overlap of colonies and patrilines on this representation was mostly due to the projection of the 34-dimensions shape space on the plane of the 2 first PCs: canonical variate analyses with cross-validation showed that the overlap between patrilines was much lower when taking more dimensions into account. Squares: 1 st patriline of each colony; circles: 2 nd patrilines; triangles: 3 rd patrilines; diamonds: 4 th patrilines. from different patrilines were significantly different (Colony 1, T² = 25.545, df num = 34, df den = 33, P = 0; Colony 2, T² = 17.108, df num = 34, df den = 33, P = 0). Correction for allometry was not required for Colony 1 as the wing size difference between the main patrilines was not significant. Allometry-corrected data confirmed that the difference in wing shape between the two main patrilines of Colony 2 was not due to their size difference (Colony 2, T² = 18.324, df num = 34, df den = 33, P = 0). These differences between colonies and between main patrilines reflected completely different shape changes (Fig 3): the angles between the vectors of these differences were not significantly similar when compared to 10,000 random vectors (θ ΔCΔP1 : 73.69°, P = 0.2096; θ ΔCΔP2 : 74.37°, P = 0.2195; θ ΔP1ΔP2 : 77.76°, P = 0.2710).

Discussion
Our results clearly show an influence of paternity-and thus of genotype-on the wing morphology of yellowjackets. In both colonies, genetic analysis revealed four different patrilines, which is congruent with previous studies on paternity in other species of Vespula yellowjackets [40] and typical of V. flavopilosa [23]. We found that the different patrilines presented small but significant differences in wing sizes. Whether this size difference is due to genetic effects or to different rearing conditions cannot be assessed yet. The nest environment of social insects, including social wasps, is remarkably homeostatic [41,42]. However, we know that the average size of cells, and thus of workers, tends to increase during the development of the colony [43]. The first worker cells built in the colony are smaller than the cells built later in the season, and the use of small cells diminishes with time in favor of larger cells. It would be interesting to assess whether patriline wing size variation may be linked to a change in sperm use patterns through time: some patrilines reared early would develop mostly in small cells while others would benefit from larger cells later. Against this environmental explanation, two studies of sperm use in yellowjackets suggest that patriline representation is constant over the course of colony development [44,45]. Regardless, the patriline influence on wing size was minimal compared to the size variation found between full-sib individuals. This confirms the results of Kovacs and collaborators [21] who found that the variance explained by patrilines was low and often nonsignificant in measurements on V. maculifrons workers. As workers only indirectly affect the colony fitness, worker size may be under less rigidly controlled genetic influence than gyne size [22].
However, our results demonstrated a clear influence of paternity on workers' wing shapes. Factors such as diet or temperature may affect the wing shape [14,15] and could vary slightly within the nest. Nevertheless, there is no evidence that individuals from different patrilines experienced different developmental conditions. These conditions would likely have an influence on the size of individuals as well [14], but the two main patrilines of Colony 1 presented no size difference. Their difference in wing shape was greater than the one between the two main patrilines of Colony 2, which differed in size as well as shape. Furthermore, even in this second colony, the differences in size between patrilines could only explain a small fraction of these shape differences. These results suggest that the father's genotype does affect the offspring's wing shapes, even if environmental factors may have contributed to the wing variation between patrilines. Unfortunately, our sample of two colonies did not allow us to calculate formal population measures such as heritability to quantify the proportions of variation due to genetic and environmental effects. On the other hand, the distinction of the two main patrilines on the first PCs within each colony suggested that genetic difference is the main cause of wing shape variation, at least within these two colonies.
The development of the wing veins has been extensively studied in Drosophila (e.g. [46][47][48]). These studies found that morphogens such as the Decapentaplegic protein (Dpp) induce the differentiation of the wing imaginal tissue into sclerotized veins. In Drosophila, the venation pattern was thought to be defined by the regulated activity of genes coding for these morphogens, such as dpp [47]. However, a recent study in Sawflies revealed that pre-patterning components directing Dpp transport were more important in defining the vein locations than the local activity of dpp [49]. How this pre-pattern is defined and under which control is still unknown [19]. Nonetheless, our results highlighted that different genotypes from a same population produced structured differences in the venation pattern: the wing venation pattern seemed under a clear genetic influence. The different patrilines were also well discriminated despite the approximately equal magnitude of shape variation between patrilines and within patrilines. This discrimination suggests that the variation within patrilines, likely induced by external factors, occurs in different directions of the shape space than the variation between patrilines, reflecting genetic effects. Finally, these differences were not equivalent, depending on the compared patrilines. Different developmental pathways may have been involved in producing the shape variation between patrilines, suggesting that the observed shape variation likely resulted from mutations across a set of genes involved in regulatory processes of the spatial activity of morphogens like Dpp [19].
Wing venation is often used for insect phylogenies based on morphology (e.g. [50]), but the venation shape assessed by geometric morphometrics is currently rarely used in this context. This could be due to methodological difficulties [51,52], though new methods are being developed to face this challenge [53,54]. Our results showed a direct link between the wing venation shape and the genotype. This link was detected despite a low genetic variation, suggesting large genetic effects. A strong genetic influence, in combination with variation in multiple independent directions, which lowers potential convergences [55,3], highlight the potential of wing shape as a phylogenetic character. Although the venation shape variation cannot be entirely attributed to the genotype [52,56], its use in phylogenetic analyses could provide additional data for species for which no genetic data are available, including fossil insects [5].
The insect wing is subject to various functional constraints concerning its size, shape and rigidity [57][58][59][60]. Consequently, the venation supporting the wing membrane likely undergoes selection on its ability to protect the wing and to control its flexibility during flight [6,55,61]. However, our results show that even small genetic differences within a population induce structured variation of the wing venation pattern. This high sensitivity of the phenotype to genetic differences suggests that the shape of wing venation may not be highly constrained by selection on the wing function: small changes in venation may not alter the flight performances of wasp workers. Further studies would be required to identify the actual influence of the wing venation shape on flight performance, and to compare the evolution of this phenotype to other features of the wing in order to detect whether they evolve under different selection regimes.
Supporting Information S1 File. Classification tables of the specimens on the basis of their wing shape using Canonical Variates Analyses with leave-one-out cross-validation. The 220 specimens of Vespula flavopilosa studied came from two colonies with 4 patrilines each. Each table presents the results when using an increasing number of Principal Components to describe the wing shape (from 2 to 34). (DOC)