Selection of the sex‐linked inhibitor of apoptosis in mountain pine beetle (Dendroctonus ponderosae) driven by enhanced expression during early overwintering

Abstract The mountain pine beetle (Dendroctonus ponderosae) is an insect native to western North America; however, its geographical range has recently expanded north in BC and east into Alberta. To understand the population structure in the areas of expansion, 16 gene‐linked microsatellites were screened and compared to neutral microsatellites using outlier analyses of F st and F ct values. One sex‐linked gene, inhibitor of apoptosis (IAP), showed a strong signature of positive selection for neo‐X alleles and was analyzed for evidence of adaptive variation. Alleles of IAP were sequenced, and differences between the neo‐X and neo‐Y alleles were consistent with neutral evolution suggesting that the neo‐Y allele may not be under functional constraints. Neo‐Y alleles were amplified from gDNA, but not effectively from cDNA, suggesting that there was little IAP expression from neo‐Y alleles. There were no differences in overall IAP expression between males and females with the common northern neo‐X allele suggesting that the neo‐X allele in males compensates for the reduced expression of neo‐Y alleles. However, males lacking the most common northern neo‐X allele thought to be selected for in northern populations had reduced overall IAP expression in early October—at a time when beetles are preparing for overwintering. This suggests that the most common allele may have more rapid upregulation. The reduced function of neo‐Y alleles of IAP suggested by both sequence differences and lower levels of expression may foster a highly selective environment for neo‐X alleles such as the common northern allele with more efficient upregulation.

−40°C are fatal to the mountain pine beetle (Carroll, Taylor, Regniere, & Safranyik, 2003). As epidemics can have severe ecological and economic impacts, it is important to model and understand mountain pine beetle outbreaks to inform future forest management practices.
The devastation caused during an epidemic is no longer limited to the historic range of D. ponderosae, but is now also a concern in new habitats where the beetles have spread in recent outbreaks. This includes movement into higher elevations in the United States, north and east in Canada, and even into a novel host, that is, jack pine (Pinus banksiana), in northern Alberta (Cullingham et al., 2011). As the mountain pine beetle has moved into these new areas, surveying the genetic diversity of beetles throughout their range has been useful in identifying both the spatial genetic structure (Samarasekera et al., 2012) and the genes potentially under selection in the expanding range (Janes et al., 2014). The identification of genes that are under selection at the northeastern limits of the beetles range may provide insights into how they have adapted to overcome the defenses of a novel host as well as the cold climates of these environments.
In looking for adaptive variation, microsatellites or Simple Tandem Repeats (STRs) are often overlooked in favor of other types of variation such as Single Nucleotide Polymorphisms (SNP's). However, microsatellites can be a useful tool to identify genes under selection that may provide local adaptation. This is especially true of Expressed Sequence Tags (ESTs) which have variable repeat lengths within the transcript which may directly result in functional differences and are easy to track through genotyping. Comparing variation in gene-linked microsatellites to the variation in neutral microsatellites within the same populations can effectively identify outliers. These outliers may be under selection in populations where certain alleles of the associated gene provide a local adaptation causing allele frequencies to be different than would be expected based on the distribution of neutral markers (Holderegger et al., 2008;Meier, Hansen, Bekkevold, Skaala, & Mensberg, 2011).
In addition to being useful markers to identify genes under selection, microsatellites may provide local adaptation by creating variation in the expression or function in different alleles of the gene (Li, Korol, Fahima, & Nevo, 2004). The expression of a gene can be greatly affected by the number of repeat units in a microsatellite (Gemayel, Cho, Boeynaems, & Verstrepen, 2012;Gymrek et al., 2016). This is particularly true if the microsatellite is in a promoter region or transcription factor binding site. Furthermore, if the microsatellite is in between two transcription factor binding sites, variation in the repeat number may interrupt or promote interactions between transcription factors (Li et al., 2004). Finally, microsatellites have been shown to inhibit nucleosome formation resulting in an open chromatin structure around microsatellites thereby promoting transcription (Gemayel et al., 2012). Variation in repeat number of a microsatellite in the coding sequence may also affect the proper formation or function of the resulting protein. In particular, microsatellites can introduce frameshift mutations, amino acid repeats can result in protein aggregation, and variable numbers of amino acids can affect protein binding interactions (Gemayel et al., 2012).
In a previous study, Samarasekera, Keeling, Bohlmann, and Murray (2011) identified 50 polymorphic gene-linked microsatellites as a resource to investigate local adaptations in mountain pine beetle populations. These were used to identify outliers showing signatures of selection. If these genes do provide a selective advantage, the location and repeat motif of the microsatellite may be able to inform what type of advantage is being conferred. As the population is expanding northeastwardly into ranges with colder climates, we hypothesize that alleles under positive selection may facilitate cold tolerance. Our screen identified the sex-linked gene for inhibitor of apoptosis (IAP) which has a microsatellite in the coding sequence as a gene showing signatures of positive selection. This gene was previously found to be upregulated overall during overwintering by mountain pine beetle larvae (Robert et al., 2016). Subsequent analysis of variation in expression levels among beetles with different genotypes during early overwintering in a northern population suggests that there are elevated levels of transcription in early overwintering in larvae with certain genotypes which may be driving this positive selection. Furthermore, the malespecific alleles of this sex-linked gene were not detected in cDNA and the sequences were divergent from neo-X linked alleles in a manner consistent with neutral evolution. Mountain pine beetle possesses a relatively recently evolved neo-X/Y sex chromosome shared with only one other species, the Jeffrey pine beetle, in the genus Dendroctonus (Zúñiga, Cisneros, Hayes, & Macias-Samano, 2002). Our findings suggest that the male-specific allele was no longer functional. Therefore, having one functional allele may have resulted in a stronger selective pressure on the neo-X alleles in male beetles. This is one example of local adaptation that may be facilitating the expansion of the range of mountain pine beetle to colder climates.

| Developing and screening a gene-linked microsatellite database in western Canada
Individual beetles collected from six sampling locations representing at least two subpopulations of MPB in Western Canada (North: Houston, BC; Mackenzie, BC; Grande Prairie, AB, and South: AB; Green Lake, BC; Banff, AB; Nancy Greene, BC [Samarasekera et al., 2012]) were genotyped at 16 gene-linked microsatellites (Table 1). Each of these microsatellites is associated with a unique scaffold (Supporting Information Table S1), and there is no evidence of linkage disequilibrium in either the neutral microsatellites (Samarasekera et al., 2012) or the gene-linked microsatellites used (Samarasekera et al., 2011). M13 tailed primers were used to amplify polymorphic microsatellites with VIC, NED, PET, or FAM fluorescent labels on the amplicons (Life Technologies Inc., Burlington, ON, Canada). The primer sequences were those developed by Samarasekera et al. (2011) with the appropriate 5′ tails. The loci were amplified using either the touchdown PCR protocol described by Samarasekera et al. (2011)  for 90 s, and 72°C for 60 s, and a final extension step at 60°C for 30 min. Fragment analysis was performed on PCR products using ABI 3130xL and GeneScan 500 LIZ size standard following the manufacturer's protocol (Life Technologies). Amplicon size was scored using GeneMapper v4.0 (Life Technologies). Genotypes were added to a database of 14 neutral microsatellite markers (Samarasekera et al., 2012) and used to determine allele frequencies, F st values, F ct values, and to identify outliers as candidates for selection using GenAlEx v6.5 (Peakall & Smouse, 2006, Arlequin v3.5 (Excoffier & Lischer, 2010), and Bayescan v2.1 (Foll & Gaggiotti, 2008). A neutral microsatellite previously identified as being sex-linked, Dpo486, was used for comparison to identify sex-linked microsatellites (Davis et al., 2009).

| Analysis of North American database
A database of 14 neutral microsatellites and one gene-linked microsatellite, MPBC6_675 (inhibitor of apoptosis, IAP), was developed using beetles from 52 populations across western North America (Boone TA B L E 1 A description of each gene-linked microsatellite including the repeat motif, location within the transcript, and the PCR method used for each of the gene-linked microsatellites screened for signatures of selection. The gene with the closest alignment from BLASTx and the organism the aligned gene is from are listed and Murray, unpublished data). This database was then divided into a female-only and male-only database because there were distinct neo-X and neo-Y linked microsatellites. The female-only database was examined for F st and F ct outliers and to determine the distribution of neo-X allele frequencies in North America. The male database was used to determine the neo-Y allele frequencies in the same 52 populations. Analysis of the genotypic data was conducted using GenAlEx v6.5 (Peakall & Smouse, 2006, Arlequin v3.5 (Excoffier & Lischer, 2010), and BayeScan v2.1 (Foll & Gaggiotti, 2008). Selected plasmid inserts were sequenced using an ABI 3130xL with T7 and SP6 primers (Life Technologies). Sequences were aligned using CodonCode Aligner (v4.2.4), and the dN/dS ratios between the sequences of the neo-X and neo-Y alleles were determined using MEGA6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013).

| Sample collection for IAP expression study
Larvae collected from Tête Juane Cache, BC in the autumn of 2008 as part of another study (Fraser, Bonnett, Keeling, & Huber, 2017) were used for expression studies of IAP. From this collection, we used larvae collected on October 3 and October 10 as these represent the time-points immediately before and after the temperature dropped below 0°C. Briefly, upon removal from the tree, larvae were immediately frozen in 1.5 ml tubes using liquid nitrogen, transported to UNBC on dry ice, and stored at −80°C until use. Samples used solely for comparison between neo-X and neo-Y expression were collected from four populations: Utah, California, Cypress Hills, and Arizona, and sent to UNBC live. Utah samples were larvae collected in June 2012, California samples were larvae collected in July 2012, Cypress Hills samples were larvae collected in November 2010, and Arizona samples were larvae, pupae, and adults collected in June 2012.

| Total RNA and DNA extraction and cDNA synthesis
The RNA and DNA from individual larvae were extracted using a

| Genotyping of IAP and prefoldin subunit 5
The microsatellite regions of the IAP and another sex-specific microsatellite region in prefoldin subunit 5 were amplified to determine the genotype of each larva and confirm its sex. These two loci were amplified in a duplex PCR using Qiagen Multiplex PCR kit (Qiagen). The final concentrations of reagents in a total volume of 10 μl were as follows: 1X for 60 s, and a final extension step at 60°C for 30 min.
Fragment analysis was performed on each of the PCR products using ABI 3130xL with GeneScan 500 LIZ size standard, scored using GeneMapper v4.0 (Life Technologies) and visually inspected to determine any miscalled alleles. Many cDNA samples were manually scored because differences in peak heights between alleles resulted in miscalled alleles. The electropherograms produced by gDNA and cDNA from each individual were compared. The ratios of the peak heights representing the neo-X and neo-Y alleles of IAP were compared between amplicons produced from gDNA and cDNA.

| Relative expression of IAP in different genotypes using RT-qPCR
The cDNA synthesized from Tête Juane Cache larvae was used as the template for RT-qPCR using an ABI 7300 system.  (Table 2).
For each target and reference gene, a standard curve was produced using three technical replicates of serially diluted pooled cDNA and used to optimize the conditions to obtain high R 2 values, efficiencies between 90 and 100%, and to determine concentrations of cDNA template within the linear dynamic range as recommended by Bustin et al. (2009). Therefore, the thermocycling conditions were different for some of the genes used; however, they all had 45 cycles of denaturation and annealing/extension (Table 3).
Individuals were divided according to genotype into four groups: females homozygous for 162, females heterozygous with a 162 allele, males with a 162 neo-X allele, and males lacking a 162 neo-X allele.
The data obtained from RT-qPCR were analyzed using REST2009 software (Qiagen; Pfaffl, Horgan, & Dempfle, 2002) to determine the relative expression of IAP between different genotypes and also between males and females.

| Sex linkage of gene-linked microsatellites
In the six Western Canadian populations genotyped, sex linkage was noted for three of the 16 gene-linked loci when compared to a previously identified genetic marker for sex linkage (Davis et al., 2009).
Two of these loci, MPBC6_675 (IAP) and MPBC6_7245 (prefoldin Step analysis. The sex linkage of IAP was consistent in beetles from 52 North American populations genotyped.

| Signatures of selection
The Western Canada microsatellite database with six sampling locations (Supporting Information Table S1)

| Analysis of spatial distribution of IAP alleles
The

| Neo-X and neo-Y allele sequences
Among the commonly observed neo-X alleles of IAP sequenced (156, 159, and 162), there was no amino acid variation outside of the microsatellite region (Figure 4). Variation was greater among the neo-Y alleles and between the neo-X and neo-Y alleles. In particular, there were two indels in the alignment of the neo-X and neo-Y sequences.
The first indel region is associated with the trinucleotide (serine) repeat and lies in a region between two predicted BIR domains. A second indel region lies between the predicted BIR2 and RING domains.
The dN/dS ratio between the 162-bp neo-X and 192-bp neo-Y alleles was 0.515 and a codon-based Z-test failed to reject neutral evolution (z = −1.44, p = 0.152). Four of the nonsynonymous changes resulting in amino acid changes were in the BIR2 and RING predicted domains.

| Genotyping of individuals
Prior to analysis of expression, individual beetles were genotyped from both gDNA and cDNA using the sex-linked microsatellites in IAP and prefoldin subunit 5 and were assigned into four groups based on their genotype. They were genotyped from both gDNA and cDNA to verify results because in some males from Utah, California, and Arizona the peak height of the neo-Y allele of IAP was not discernible above background when genotyped from cDNA alone due to a lack of expression.

F I G U R E 2
The F st and F ct values for neutral and gene-linked microsatellites genotyped in mountain pine beetles collected in Western Canada. * indicates an outlier based on both F st and F ct values (p < 0.01)

| Qualitative relative expression of IAP alleles
The electropherograms from the amplification of prefoldin subunit 5 produced peaks with equal neo-X to neo-Y peak height ratios from gDNA and cDNA. In contrast, the electropherograms from the amplification of IAP from cDNA produced peaks with higher neo-X to neo-Y peak height ratios than those from gDNA ( Figure 5). Based on the electropherograms from 18 male samples from Tête Juane Cache, BC the mean ratio of neo-X to neo-Y allele peak heights from cDNA was an average of 9.55 times higher than the ratio from gDNA (p < 0.005). A similar trend was found in all male individuals from Utah, California, Arizona, and Cypress Hills. In each location, the peak produced by the neo-Y allele was reduced in the cDNA compared to the gDNA. However, the ratio was not quantifiable since the peak corresponding to the neo-Y allele was not visible above background in many of the electropherograms.

| Quantitative relative expression of IAP in larvae
The optimized conditions for RT-qPCR using cDNA from Tête Juane isolated SNP haplogroups (Bracewell, Bentz, Sullivan, & Good, 2017;Dowle et al., 2017). Using this single gene, we find consistency with the populations from the western United States being similar; however, in Canada, eastern US, and central US, we find more F I G U R E 3 Spatial distribution of (a) neo-X alleles and (b) neo-Y alleles of IAP genotyped in mountain pine beetles collected from the indicated sampling locations across western North America differences among neo-Y alleles within populations than would be predicted by geographic range of the SNP haplogroups (Figure 3b).

This suggests that the inclusion of neo-Y microsatellite markers with
SNPs Dowle et al., 2017) will provide a finer resolution of within haplogroup variation which can inform evolution of the neo-Y chromosome in these areas.
The independent evolution of neo-X and neo-Y alleles in D. ponderosae is of particular interest because there has been a recent fusion between an ancestral autosomal and X chromosome as well as the loss of the ancestral Y chromosome (Keeling et al., 2013;Lanier & Wood, 1968

F I G U R E 5
The peak heights of neo-X and neo-Y alleles of both IAP and prefoldin subunit 5 amplified from gDNA compared to the same alleles amplified from cDNA from a single individual larva. The neo-Y allele of IAP amplified from cDNA is greatly reduced compared to the allele amplified from gDNA, but no difference is observable for the neo-Y allele of prefoldin subunit 5 two sex chromosomes (X and Y), whereas D. ponderosae and D. jeffreyi have 11 autosomal chromosomes and two sex chromosomes (neo-X and neo-Y) (Zúñiga et al., 2002). It is expected that over time sequence divergence will occur between the ancestral autosomal portions of the neo-X and neo-Y chromosomes explaining the distinct IAP alleles observed. However, it is still early in this divergence so we were able to amplify and study both the neo-X and neo-Y alleles with common primer sets.

Analysis of sequence variation between the neo-X and neo-Y
IAP alleles failed to reject neutral evolution. This variation suggests that the neo-Y IAP may no longer be evolving under functional constraints indicating that it may not contribute substantially to overall IAP production. Furthermore, when male samples were genotyped the electropherograms for IAP amplified from cDNA had a reduced peak height produced by the neo-Y allele. The peak height ratio of the IAP neo-X allele compared to the neo-Y allele was consistently larger in the electropherograms produced from amplification of cDNA than those from gDNA across all populations, representing seven different neo-Y alleles ( Figure 5). As the ratio of neo-X to neo-Y alleles in genomic DNA is one to one, this suggests that there were more neo-X allele transcripts of IAP than neo-Y allele transcripts to acts as templates for cDNA synthesis due to reduced expression of the neo-Y allele. Reduced expression may have led to the loss of functional constraints which resulted in sequence variation between the neo-Y and neo-X alleles being consistent with neutral evolution.
In other insects, such as D. miranda with neo-X/neo-Y sex chromosomes, the neo-Y chromosome is degenerating (Bachtrog & Charlesworth, 2002). Based on the sequence variation between the IAP neo-X and neo-Y alleles presented in this study, there is likely degeneration in the D. ponderosae neo-Y chromosome. This is supported by recent findings that there has been degeneration in the neo-Y chromosome of D. ponderosae. Bracewell et al. (2017) report large deletions in the neo-Y chromosome suggesting degeneration, whereas we find divergence usually associated with a pseudogene (equal number of synonymous and nonsynonymous changes) within a single gene between the neo-X and neo-Y chromosomes. These pieces of evidence point to neo-Y degeneration which may provide a highly selective environment for the orthologous neo-X alleles.
Outside of the variable number of serine repeats, there was no sequence variation in the coding sequences of neo-X alleles analyzed, suggesting that differences in protein activity of different neo-X alleles are unlikely. Therefore, potential expression differences among neo-X alleles were examined.
Overall IAP expression among all larvae was found to increase significantly from 3 October to 10 October consistent with an overall trend of increasing IAP expression in larvae preparing for overwintering which had a 3.13-fold increase from 26 September to 7 November 2008 in the same collection (Robert et al., 2016). This increase in expression as the larvae prepare for overwintering suggests that IAP may play a role in cold tolerance. This potential role in cold tolerance is supported by the role of IAP as a caspase inhibitor as apoptosis involving caspases has been shown to be inducible by cold shock (Fransen, Dieker, Hildebrands, Berden, & van der Vlag, 2011). Furthermore, in Drosophila melanogaster, an inhibitor of apoptosis Bcl-2 was shown to be upregulated at the protein level during rapid cold-hardening to compensate for dysregulated apoptosis which is common during cold shock (Yi, Moore, & Lee, 2007) Therefore, elevated levels of IAP expression may help the larvae to survive overwintering.
The relative expression of IAP in overwintering larvae was com- The lack of expression differences between males and females suggests that dosage compensation may be occurring. As discussed, it appears that the neo-Y alleles have reduced expression based on the electropherograms from cDNA. Therefore, dosage compensation would be necessary to allow equal expression of IAP in males and females. Although some form of dosage compensation is known to occur in many species with both neo-X/neo-Y and X/Y sex chromosomes, it has not been reported in bark beetles. D. miranda has a neo-X/neo-Y sex chromosome system which has dosage compensation partly developed for neo-X alleles in males (Bachtrog & Charlesworth, 2002). In particular, there is increased expression of genes on the neo-X chromosome which have degenerated orthologous genes on the neo-Y chromosome (Marín, Siegal, & Baker, 2000). In D. ponderosae, the sequence variation and reduced expression suggests that the neo-Y IAP allele may be degenerating.
Furthermore, another sex-linked gene, prefoldin subunit 5, did not show variation in the peak heights of neo-X and neo-Y alleles. This is also consistent with the neo-X/neo-Y dosage compensation in D. miranda which is not chromosome wide (Bone & Kuroda, 1996).
Although there is little known about dosage compensation in Coleoptera, Prince, Kirkland, and Demuth (2010) found that the majority of X-linked genes in the red flour beetle, Tribolium castaneum, have increased expression in females. As the X-linked genes in male flour beetles were expressed at approximately the same levels as autosomal genes, the increased expression of X-linked genes in females is believed to be an evolutionary side effect of an imperfect dosage compensation mechanism. This hyperexpression of X-linked genes shows that differential expression of sex chromosomes occurs in Coleoptera. However, the equal expression of IAP in male and female D. ponderosae suggests that, in contrast to T. castaneum, it is unlikely that neo-X linked genes are globally hyperexpressed in females. It is possible that D. ponderosae is more similar to Xenos vesparum in which the chromosomal segment corresponding to the ancestral X shows dosage compensation but dosage compensation in the recently added segment of the X chromosome is incomplete (Mahajan & Bachtrog, 2015). However, the expression of more neo-X linked genes particularly those in the ancestral X segment would need to be determined to make a definite conclusion on the extent and mode of dosage compensation in D. ponderosae.
Despite recent studies on sex chromosome evolution in D. ponderosae Dowle et al., 2017), there is still little known about dosage compensation and much to be explored on a finer scale regarding sex chromosome evolution in bark beetles, so studying variation in both the sequence and expression of neo-X and neo-Y linked genes can help determine how the sex chromosomes are evolving and predict how they will continue to evolve. The sequence variation between the IAP neo-X and neo-Y alleles suggests that the neo-Y chromosome is degenerating in D. ponderosae. Such degeneration fosters a highly selective environment for neo-X alleles as they must compensate for an orthologous neo-Y allele with a reduced or lost function. In the case of IAP, the 162-bp allele appears to have faster upregulation compared to the other neo-X alleles, which could provide a selective advantage allowing mountain pine beetle larvae in the north to survive early cold snaps.
The sequence variation and expression differences in mountain pine beetle IAP suggest that the neo-Y allele is not contributing substantially to the overall production of the protein. Experience Awards to LCH and GKK.

CO N FLI C T O F I NTE R E S T
None declared.