\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Chromosome-level genome assembly reveals genomic architecture of
northern range expansion in the mountain pine beetle, \emph{Dendroctonus
ponderosae} Hopkins (Coleoptera: Curculionidae)}



\author[1]{Christopher I. Keeling}%
\author[2]{Erin O. Campbell}%
\author[3]{Philip D. Batista}%
\author[2]{Victor A. Shegelski}%
\author[2]{Stephen A. L. Trevoy}%
\author[3]{Dezene P. W. Huber}%
\author[4]{Jasmine K. Janes}%
\author[2]{Felix A. H. Sperling}%
\affil[1]{Laurentian Forestry Centre}%
\affil[2]{University of Alberta}%
\affil[3]{University of Northern British Columbia}%
\affil[4]{Vancouver Island University}%


\vspace{-1em}



  
  \date{May 25, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Genome sequencing methods and assembly tools have improved dramatically
since the 2013 publication of draft genome assemblies for the mountain
pine beetle, \emph{Dendroctonus ponderosae} Hopkins (Coleoptera:
Curculionidae). We conducted proximity ligation library sequencing and
scaffolding to improve contiguity, and then used linkage mapping and
recent bioinformatic tools for correction and further improvement. The
new assemblies have dramatically improved contiguity and gaps compared
to the originals: N50 values increased 26- to 36-fold, and the number of
gaps were reduced by half. Ninety percent of the content of the
assemblies is now contained in 12 and 11 scaffolds for the female and
male assemblies, respectively. Based on linkage mapping information, the
12 largest scaffolds in both assemblies represent all 11 autosomal
chromosomes and the neo-X chromosome. These assemblies now have nearly
chromosome-sized scaffolds and will be instrumental for studying genomic
architecture, chromosome evolution, population genomics, functional
genomics, and adaptation in this and other pest insects. We also
identified regions in two chromosomes, including the ancestral-X portion
of the neo-X chromosome, with elevated differentiation between northern
and southern Canadian populations.%
\end{abstract}%



\sloppy


\textbf{Chromosome-level genome assembly reveals genomic architecture of
northern range expansion in the mountain pine beetle, \emph{Dendroctonus
ponderosae} Hopkins (Coleoptera: Curculionidae)}

Christopher I. Keeling\textsuperscript{1,2,*}, Erin O.
Campbell\textsuperscript{3}, Philip D. Batista\textsuperscript{4},
Victor A. Shegelski\textsuperscript{3}, Stephen A. L.
Trevoy\textsuperscript{3}, Dezene P. W. Huber\textsuperscript{4},
Jasmine K. Janes\textsuperscript{5,6}, and Felix A. H.
Sperling\textsuperscript{3,*}

1. Laurentian Forestry Centre, Canadian Forest Service, Natural
Resources Canada, Qu\selectlanguage{ngerman}ébec, QC, Canada

2. Département de biochimie, de microbiologie et de bio-informatique,
Université Laval, Québec, QC, Canada

3. Department of Biological Sciences, University of Alberta, Edmonton,
AB, Canada

4. Faculty of Environment, University of Northern British Columbia,
Prince George, BC, Canada

5. Biology Department, Vancouver Island University, Nanaimo, BC, Canada

6. School of Environmental and Rural Studies, University of New England,
Armidale, NSW, Australia

Running title (45 characters max): Genome assemblies of the mountain
pine beetle

ORCID iDs:

CIK, 0000-0001-7885-8558; EOC, 0000-0001-8546-0636; PDB,
0000-0002-9598-2793; VAS, 0000-0003-2438-9447; SALT,
0000-0003-1217-406X; DPWH, 0000-0002-6495-1759; JKJ,
0000-0002-4511-2087; FAHS, 0000-0001-5148-4226.

* Corresponding authors: christopher.keeling@canada.ca;
felix.sperling@ualberta.ca



\textbf{Abstract}

Genome sequencing methods and assembly tools have improved dramatically
since the 2013 publication of draft genome assemblies for the mountain
pine beetle, \emph{Dendroctonus ponderosae} Hopkins (Coleoptera:
Curculionidae). We conducted proximity ligation library sequencing and
scaffolding to improve contiguity, and then used linkage mapping and
recent bioinformatic tools for correction and further improvement. The
new assemblies have dramatically improved contiguity and gaps compared
to the originals: N50 values increased 26- to 36-fold, and the number of
gaps were reduced by half. Ninety percent of the content of the
assemblies is now contained in 12 and 11 scaffolds for the female and
male assemblies, respectively. Based on linkage mapping information, the
12 largest scaffolds in both assemblies represent all 11 autosomal
chromosomes and the neo-X chromosome. These assemblies now have nearly
chromosome-sized scaffolds and will be instrumental for studying genomic
architecture, chromosome evolution, population genomics, functional
genomics, and adaptation in this and other pest insects. We also
identified regions in two chromosomes, including the ancestral-X portion
of the neo-X chromosome, with elevated differentiation between northern
and southern Canadian populations.

\textbf{Keywords}

proximity ligation sequencing, Hi-C, linkage map, chromosome
architecture, sex chromosomes, forest pest management, insect ecology

\textbf{Introduction}

The exponential increase in the availability of whole genomes over the
past decade has largely been derived from short-read sequences. However,
while short reads can be assembled into larger contigs and scaffolds,
sequencing gaps and repetitive regions often prevent full contiguous
assembly and correct ordering (Fierst, 2015). Many genomes assembled
from short reads remain highly fragmented and are best considered draft
versions. Long-read sequencing technologies can overcome some issues
associated with gaps and repetitive regions, but they suffer from higher
sequencing error rates (Bleidorn, 2015; Rhoads \& Au, 2015). Chromosome
conformation capture sequencing, or Hi-C, is a method to study
three-dimensional folding of chromosomes that offers an elegant solution
to these challenges. Specifically, the Hi-C method detects closely
linked sections (approximately 10-100 Kb) of DNA via chromatin contacts
(van Berkum et al., 2010), producing highly accurate information on
chromosome-scale genomic position (Cairns et al., 2016). When the Hi-C
approach is coupled with an existing short-read draft genome,
scaffolding is dramatically improved by additional anchor points
interspersed throughout the short-read assembly (Li et al., 2018).

The availability of high-quality whole genomes facilitates an enhanced
understanding of genomic organization, the mechanisms that lead to
deviations from `common' architectural patterns, and the potential
outcomes of such deviations. For example, a new, highly-contiguous
genome assembly for the honey bee, \emph{Apis mellifera} L.
(Hymenoptera: Apidae), supported the detection of inversions associated
with local adaptation in high altitude populations (Christmas et al.,
2019). Comparative genomics has also contributed to a growing
recognition of genomic islands of divergence and speciation (Gagnaire et
al., 2013; Renaut et al., 2013). Numerous recent studies highlight the
ecological roles that divergent regions within genomes may have and how
they are formed (e.g., Bay \& Ruegg, 2017; Ma et al., 2018; Renaut et
al., 2013; Tavares et al., 2018).

Neo-sex chromosomes and the degeneration of male sex chromosomes (e.g.,
Muller's ratchet), present other aspects of genomic speciation that can
be explored using high-quality genome assemblies. For example,
assessments of neo-sex chromosome formation and prevalence, and the
subsequent deterioration of male sex chromosomes, can provide insight
into the mechanisms and consequences of their evolution. High-quality
insect genomes are increasingly valuable because they present numerous
complex sex-determination mechanisms that can affect their stability
through recombination, mutation, and drift (Blackmon et al., 2016).
Recent research has highlighted the interaction between sex chromosomes
and reproductive isolation (Bracewell et al., 2017), and insect genomes
with newly derived neo-sex chromosomes may be better suited for studying
male sex chromosome degeneration than mammal genomes.

Genomic resources for mountain pine beetle, \emph{Dendroctonus
ponderosae} Hopkins (Coleoptera: Curculionidae), have grown considerably
in recent years (Cullingham et al., 2018; Keeling et al., 2012; Keeling
et al., 2013c). Much of this research and data growth has been in
response to ongoing massive population outbreaks in western North
America driven by climate change, which have resulted in severe forest
disturbances and losses (e.g., Raffa et al., 2008; Saab et al., 2014;
Sambaraju \& Goodsman, 2021). Genetics and genomics were quickly
recognized as tools to help model and manage outbreaks by better
understanding population dynamics (Goodsman et al., 2016; James et al.,
2016); signatures of selection (Batista et al., 2016; Janes et al.,
2014); the processes of olfaction (Andersson et al., 2013; Andersson et
al., 2019; Chiu et al., 2019b; Keeling et al., 2013b), pheromone
biosynthesis (Aw et al., 2010; Chiu et al., 2018, 2019a; Keeling et al.,
2013a; Keeling et al., 2016; Nadeau et al., 2017), host colonization
(Chiu et al., 2019c; Robert et al., 2013) and overwintering (Robert et
al., 2016); and the coevolution between the beetle, its hosts, and its
symbiotic microorganisms (James et al., 2011).

While recent focus has been on applied pest management outcomes
for\emph{D. ponderosae} , there are numerous theoretical benefits as
well, and the two often complement each other. For example, the draft
mountain pine beetle genome (Keeling et al., 2013c) revealed
considerable synteny with \emph{Tribolium castaneum} (Herbst)
(Coleoptera: Tenebrionidae), the most closely related genome assembly at
the time. That work highlighted the important roles that specific gene
families and the neo-XY chromosomal system may have in contributing to
the ecological success of the mountain pine beetle (Keeling et al.,
2013c). That study also noted the lack of linkage map data to further
refine our knowledge of the location and formation of the neo-XY system.
Traditionally, linkage maps have been used to quantify effect sizes that
certain loci have on traits of interest (e.g., quantitative trait loci)
(Berlin et al., 2017), which can greatly accelerate breeding programs
and our understanding of evolutionary processes and outcomes. However,
linkage maps can also assist in genome assembly of non-model organisms
(Bartholomé et al., 2014) and provide insight into comparative analyses
of genomic architecture (i.e., synteny, collinearity, chromosomal
rearrangements, etc.) (Barth et al., 2019; Butler et al., 2017).

The karyotype of \emph{Dendroctonus} spp. (Coleoptera: Curculionidae)
varies from 5 autosome pairs (AA) + neo-XY to 14 AA +
Xy\textsubscript{p} (Zúñiga et al., 2002), with males being the
heterogametic sex in all cases. In \emph{D. ponderosae} and \emph{D.
jeffreyi} Hopkins, the karyotype is 11 AA + neo-XY. Neo-sex chromosomes
result from the fusion of a sex chromosome with an autosome. Three
other\emph{Dendroctonus} spp. have a neo-XY karyotype: \emph{D.
adjunctus}Blandford (6 AA + neo-XY), \emph{D. approximatus} Dietz (5 AA
+ neo-XY), and \emph{D. brevicomis} LeConte (5 AA + neo-XY). It is not
known whether the neo-XY chromosomes in these species are orthologous to
those in \emph{D. ponderosae} and \emph{D. jefferyi} . In \emph{D.
ponderosae} , the neo-X is believed to have originated from an ancestral
12 AA + Xy\textsubscript{p} karyotype through the fusion of the
ancestral-X chromosome with one copy of the largest autosome (Lanier,
1981). The ancestral-Y (y\textsubscript{p}) chromosome was lost, to be
replaced by the other copy of the fused autosome, becoming neo-Y.
However, several species in this genus have the basal 14 AA +
Xy\textsubscript{p}karyotype (Zúñiga et al., 2002), suggesting that the
large autosome that became part of the neo-sex chromosomes in \emph{D.
ponderosae} may itself have been the fusion of three ancestral
autosomes. Thus, a large portion of the genome content of \emph{D.
ponderosae} is in the sex chromosomes. Apart from understanding the
genome architecture, delineating the sex chromosomes and the genes they
contain will help us to understand how \emph{D. ponderosae} has become
one of the most important pest species in its genus, and a keystone
species in western North American pine forests.

Here we provide chromosome-level male and female genome assemblies based
on new proximity ligation and linkage mapping data to further support
research on mountain pine beetle, which continues to exert major
ecological and economic effects on pine forests in Canada and the USA.
We additionally use whole-genome pool-seq and SNP data to assess
signatures of natural selection and investigate genomic linkage in the
context of rapid population, geographic, and host range expansions, as
well as neo-sex chromosome evolution. Our results broadly characterize
the neo-sex chromosomes and indicate that local adaptation in Canadian
mountain pine beetle populations appears to be largely mediated by one
autosomal region and the ancestral-X portion of the neo-X chromosome.

\textbf{Materials and Methods}

\emph{Genome assembly and annotation}

\emph{Insects for genome scaffolding}

Female adult beetles were collected in August 2017 from a lodgepole pine
forest approximately 60 km south of Grande Prairie, Alberta, Canada,
54.6942 N 118.6694 W. Specimens from this collection have been vouchered
at the E. H. Strickland Entomological Museum at the University of
Alberta with UASM numbers 391986-391987.

\emph{Chicago library preparation and sequencing}

Dovetail Genomics LLC prepared a Chicago library as described previously
(Putnam et al., 2016). Briefly, approximately 500 ng of high molecular
weight gDNA (mean fragment length = 80 kb) from an individual female
beetle (NCBI BioSample: SAMN14918906) was reconstituted into
chromatin\emph{in vitro} and fixed with formaldehyde. Fixed chromatin
was digested with \emph{DpnII} , the 5'-overhangs filled in with
biotinylated nucleotides, and then free blunt ends were ligated. After
ligation, crosslinks were reversed, and the DNA was purified from
associated protein. Purified DNA was treated to remove biotin that was
not internal to ligated fragments. The DNA was then sheared to
approximately 350 bp mean fragment length and sequencing libraries were
generated using NEBNext Ultra enzymes and Illumina-compatible adapters.
Biotin-containing fragments were isolated using streptavidin beads
before PCR enrichment of the library. The library was sequenced on an
Illumina HiSeq X platform to produce 510 million 2x150 bp paired-end
reads (NCBI SRA: SRR11965410), which provided 747 x physical coverage of
the genome (1-100 kb pairs).

\emph{Hi-C library preparation and sequencing}

Dovetail Genomics LLC prepared a Hi-C library in a manner similar to
that described previously (Lieberman-Aiden et al., 2009). Briefly,
chromatin from another individual female beetle (NCBI BioSample:
SAMN14918916) was fixed in place with formaldehyde in the nucleus and
then extracted fixed chromatin was digested with \emph{DpnII} . The 5'
overhangs were filled in with biotinylated nucleotides and then the free
blunt ends were ligated. After ligation, crosslinks were reversed, the
DNA was purified from protein, and a library was prepared as above and
sequenced on an Illumina HiSeq X platform to produce 514 million 2x150
bp paired-end reads (NCBI SRA: SRR11965415), which provided 37\selectlanguage{english} 483 x
physical coverage of the genome (10-10\selectlanguage{english} 000 kb pairs).

\emph{Scaffolding the assembly with HiRise}

A schematic for the assembly process is shown in Fig. 1. We used the
pre-existing male (NCBI BioProject: PRJNA162621) and female (NCBI
BioProject: PRJNA179493) genome assemblies (Keeling et al., 2013c) as
the input \emph{de novo} assemblies for scaffolding, along with their
underlying shotgun reads (NCBI SRA: SRR546181 and SRR546185 for male
assembly, and SRR546193 for female assembly).

Separately for each sex, the \emph{de novo} assembly, shotgun reads,
Chicago library reads, and Dovetail Hi-C library reads were used as
input data for HiRise, Dovetail's proprietary software pipeline designed
specifically for using proximity ligation data to scaffold genome
assemblies (Putnam et al., 2016). First, shotgun and Chicago library
sequences were aligned to the draft input assembly using the SNAP
v1.0beta read mapper (Zaharia et al., 2011). Separations of Chicago read
pairs mapped within draft scaffolds were analyzed by HiRise to produce a
likelihood model for genomic distance between read pairs; the model was
used to identify and break putative mis-joins, to score prospective
joins, and make joins. Second, Hi-C library sequences were aligned and
scaffolded following the same method. Third, shotgun sequences were used
to close gaps between contigs.

\emph{Correction and Improvement of HiRise assemblies}

We compared the male and female HiRise assemblies to each other using
NUCmer (part of MUMmer v 3.23\_3, Kurtz et al., 2004), and to linkage
map information (described below) using ALLMAPS v0.8.12 (Tang et al.,
2015), to identify inconsistencies between the two assemblies, and
between the assemblies and the linkage map information. Aided by
NUCmer-based mummerplot comparisons of the relevant scaffolds, we
resolved inconsistencies by manually flipping and/or shuffling the
Chicago scaffolds that appeared to be mis-joined in the final HiRise
step such that both assemblies were mutually consistent and were also
consistent with the linkage map information. When there were
inconsistencies between the male and female assemblies, we deemed the
assembly most consistent with the linkage map information to be correct.
At this stage we also shortened gaps greater than 1000 bp that were in
the original draft assembly, to 999 bp. We then scaffolded and
gap-filled the assemblies with the corresponding assembly from the other
sex using LINKS v1.8.7 (Warren et al., 2015) and RAILS v1.4.1 (Warren,
2016). We filled the remaining gaps in the scaffolds with ABySS-Sealer
v2.1.5\_1 (Paulino et al., 2015) using shotgun sequence data from the
original \emph{de novo} assembly (NCBI SRA: SRR546193 for female
assembly, and SRR546181 and SRR546185 for male assembly). We removed
contaminant scaffolds identified by PhylOligo v1.0 (Mallet et al.,
2017). We used BUSCO v4.1.4 (Simao et al., 2015) analysis with the
insecta\_odb10.2020-09-10 lineage dataset (Kriventseva et al., 2019) to
assess improvement of gene content on the final genome assemblies, and
the input \emph{de novo} assemblies. Assembly consistency plots were
generated using JupiterPlot v1.0 (Chu, 2020) with Circos v0.69-3
(Krzywinski et al., 2009).

\emph{Gene prediction and annotation of genome assemblies}

We modeled repetitive sequences with RepeatModeler v2.0 (part of Dfam TE
Tools v1.0, Dfam Consortium, 2019) and masked the assemblies of these
repetitive regions with RepeatMasker v4.0.7\_4 (Smit et al., 2020). We
then completed three rounds of Maker v2.31.10\_2 (Holt \& Yandell, 2011)
gene prediction on scaffolds longer than 5 kb. For the first round, we
used NCBI Coleopteran RefSeq proteins (downloaded 2019-11-12) for
protein homology evidence, and NCBI \emph{D. ponderosae} ESTs
(downloaded 2019-11-12); \emph{D. ponderosae} GenBank FLcDNAs:
BT126413-BT128693 and JQ855638-JQ855707; and NCBI
\emph{Dendroctonus}spp. TSA accessions: GABX00000000, GAFI00000000,
GAFW00000000, GAFX00000000, GDAR00000000, and GGKQ00000000 for EST
evidence. Alignments were completed using exonerate v2.4.0 (Slater \&
Birney, 2005) and BLAST+ 2.9.0 (Camacho et al., 2009). The results from
the first round were used with BUSCO v3.0.2\_2 to train AUGUSTUS v3.3.3
(Stanke et al., 2006). Round two used AUGUSTUS, SNAP v 2013-11-29\_1
(Korf, 2004), and GeneMark-ES v2019-05 (Lomsadze et al., 2005)
for\emph{ab initio} gene prediction. We used the results from round two
to directly train AUGUSTUS for the final round of Maker along with SNAP
and GeneMark-ES again. We then used InterProScan v5.40-77.0 (Jones et
al., 2014) and BLAST+ v2.10.0 against the UniProtKB/Swiss-Prot 2020\_01
database to functionally annotate the predicted proteins from the gene
predictions. BUSCO v4.1.4 analysis with the insecta\_odb10.2020-09-10
lineage dataset was completed on the predicted proteins from the first
transcript variant of each gene model. The annotated assemblies are
available with the following NCBI accessions: BioProjects, PRJNA638274
(male) and PRJNA638278 (female); genome accessions JAFETF000000000
(male) and JAFETG000000000 (female).

\emph{Linkage Mapping}

\emph{Controlled crosses}

We collected mountain pine beetle larvae (F0 generation) from infected
lodgepole pine trees from Smokey River Lowlands near Grande Prairie
(Alberta, Canada; 54deg21.376' N, 118deg19.112' W) and BURNCO Quarry near
Canmore (Alberta, Canada; 51deg04.026' N, 115deg17.237' W) in October of
2014 (Supp. Fig. 1). These sites represent the two genetically distinct
sub-populations previously identified within Canada (Batista et al.,
2016; Janes et al., 2014; Samarasekera et al., 2012) and are referred to
in this study as North and South, respectively. Collected larvae
continued development in their natal bolts at 4degC until spring 2015.
Next, we placed bolts in opaque emergence boxes at room temperature
(approximately 22degC) and mature mountain pine beetles were collected
daily for a maximum of 15 days. Freshly emerged mountain pine beetles
were sexed using the auditory method described by Rosenberger et al.
(2016) to minimize physical damage. We temporarily stored the emerged
adults at 4degC until a sufficient number had been collected to begin
crosses (maximum five days). Further details are provided in Trevoy et
al. (2019).

We reared North x South crosses of the F0 generation in 55 cm bolts from
three lodgepole pine, \emph{Pinus contorta} Douglas ex Louden (Pinales:
Pinaceae), trees felled in May 2015 at Nojack, Alberta (53deg36.103' N,
115deg35.239' W). Females were placed in 1.5 ml microcentrifuge tubes
attached to bolts to encourage gallery construction within a specific
bolt. Males from the opposite population were placed in the tube 48 h
later. Females or males that failed to enter the bolt within 24 h were
replaced with alternates from the same site of origin. We stored F0
crosses in a sealed room at approximately 22degC for six weeks, at which
point we transferred them to emergence boxes. F1 progeny were collected
over 15 days and physically sexed using the non-destructive auditory
method. Mountain pine beetle adults have been known to establish
secondary galleries in different trees (Janes et al., 2016). To ensure
that F1 beetles were assigned to the correct gallery, and therefore the
correct cross, each bolt was stripped and visually inspected. Any
progeny from a bolt containing multiple galleries was discarded from
further analysis.

We crossed F1 beetles with gallery siblings to produce highly related F2
individuals. Males tend to guard the entrance of a gallery and pack it
with frass (Reid, 1962). Thus, after seven days, we peeled back the wood
around each gallery to a depth of 5 cm and the male was removed for
genotyping. In accordance with Amman (1972), females were expected to
complete gallery construction within three weeks. Thus, we extracted
females at a later date in a similar fashion to males. Bolts were placed
in emergence boxes at approximately 22degC until additional F2 progeny
were no longer recovered. All F2 individuals were stored at --20degC prior
to DNA extraction after physical sexing, using the dimorphic seventh
tergite (Safranyik \& Carroll, 2006). Sexing using the dimorphic seventh
tergite was reserved for F2 individuals because, while more accurate
than the auditory method, it typically causes damage to the specimen.
The crossing design and family summaries are illustrated in Supp. Fig.
1. Specimens have been vouchered at the E. H. Strickland Entomological
Museum at the University of Alberta with UASM numbers 391988-391991.

\emph{DNA extraction, sequencing, and trimming of linkage mapping
samples}

We extracted DNA from 229 F1 and F2 individuals from the North X South
crosses using the Qiagen DNeasy Blood and Tissue kit with optional RNase
treatment. Samples underwent SNP library preparation following the ddRAD
protocol of Peterson et al. (2012) using \emph{PstI} and
\emph{MspI}restriction enzymes. Illumina 75 bp single-end sequences were
produced on a NextSeq500 at the University of Alberta Molecular Biology
Services Unit. We demultiplexed the resulting sequence reads using
the\emph{process\_radtags} module of Stacks v2.0 (Catchen et al., 2011;
Rochette et al., 2019). Remnant Illumina adapters and the \emph{PstI}cut
site at the 5' end of each read were trimmed using Cutadapt v1.10
(Martin, 2011). The trimmed length for each read was 62 bp. We have
deposited the demultiplexed, trimmed reads for these samples (NCBI
BioSamples SAMN17498181-SAMN17498396) into the NCBI Sequence Read
Archive under accession PRJNA694171.

\emph{Data processing and linkage map construction in Lep-MAP3}

We used mpileup in SAMtools v1.9 (Li et al., 2009), and the
pileupParser2 and pileup2posterior scripts implemented in Lep-MAP3
(Rastas, 2017), to align trimmed reads to both the male and female Hi-C
HiRise genome assemblies separately. Resulting posterior files were used
as input into Lep-MAP3 to produce male- and female-aligned linkage maps.
We used identity-by-descent (IBD) scores to verify the assignment of
individuals to discrete families by removing individuals with less than
25\% IBD to at least half of the individuals within their respective
families. The ParentCall2 module imputed missing genotypes in the F1
parents that were not recovered from the bolts. The Filtering2 module
removed markers with high segregation distortion or excessive missing
data (data tolerance score of 0.01, following Lep-MAP3 recommendations).

Next, we used the SeparateChromosomes2 module of Lep-MAP3 to separate
SNPs into distinct linkage groups representing putative chromosomes. We
required retained linkage groups to contain at least 70 SNPs, and set
the informativeMask parameter to 23, which excluded markers that were
informative only for the fathers (i.e., we retained markers that were
either informative for the mothers, or for both mothers and fathers).
Including markers informative only for fathers substantially reduced the
number of SNPs assigned to linkage groups. We adjusted LOD scores until
the number of retained linkage groups closely matched the known number
of chromosomes (11 autosomes + 2 neo-sex chromosomes) based on mountain
pine beetle karyology (Lanier \& Wood, 1968). In general, the
appropriate LOD score should be similar to the number of chromosomes in
the genome (Rastas, 2017).

We then ordered each linkage group in both maps five times using the
OrderMarkers2 module of Lep-MAP3 and selected the marker order for each
group with the highest likelihood score. We checked each file for
incorrect marker ordering by visualizing linkage group graphs with xdot
v1.1 (Fonseca, 2019). If any of these graphs indicated improper marker
ordering, we discarded that replicate, chose the replicate with the
next-highest likelihood score, and checked it again. This produced
separate SNP recombination distances for male and female specimens in
each linkage group, which were used as input for ALLMAPS (Tang et al.,
2015). These linkage maps were used to inform joining, ordering, and
orientation of scaffolds in the male and female Hi-C HiRise genome
assemblies (described above). After these assembly modifications and
subsequent steps were completed, we reproduced the male- and
female-aligned linkage maps and ALLMAPS figures using the final versions
of the male and female genome assemblies and the same parameters
described above. We then repeated the marker ordering step in Lep-MAP3,
this time outputting a single, sex-averaged distance for each SNP in the
linkage groups in order to produce chromosome maps using the
LinkageMapView package (Ouellette et al., 2018) in R v3.6.1 (R Core
Team, 2020). We visualized each chromosome as a density map to identify
regions with strong genetic linkage, indicated by shorter per-locus cM
distances.

\emph{Assessment of local adaptation in Canadian mountain pine beetle
populations}

\emph{Mapping}

Paired-end pool-seq reads described previously (Keeling et al., 2013c)
were obtained from the NCBI SRA database (Fairview (FV) - SRR073440;
Kananaskis (KA) - SRR086167; Terrace (TR) -SRR073431; Whitecourt (WC) -
SRR073441; Cypress Hills (CH) - SRR086168; Valhalla (VA) - SRR086169;
Houston (HO) - SRR086170). Paired-end reads were trimmed using BBDuk
v37.25 (Bushnell, 2020) (trimq = 20, minlength = 10) prior to mapping.
Two datasets containing northern (HO, TR, WC, FV) and southern (VA, KA,
CH) Canadian populations were created and mapped to the final female
genome assembly using the BBMap v36.92 (Bushnell, 2020) plugin for
Geneious v10.1.3 (Biomatters, Auckland, New Zealand) with default
settings.

\emph{Differentiation and Selection}

We used Popoolation2 (Kofler et al., 2011) to evaluate nucleotide
diversity (\selectlanguage{greek}π\selectlanguage{english}), differentiation (\emph{F\textsubscript{ST}} ) and
Tajima's\emph{D} between northern and southern pool-seq
datasets.\emph{F\textsubscript{ST}} values were calculated using a
non-overlapping sliding window approach (min-count -2, min-coverage 4,
max-coverage 70, min-covered-fraction 0.2, pool-size 30, window-size
10 000, step-size 10 000). Tajima's \emph{D} values were calculated with
Npstat v1 (Ferretti et al., 2013) using a sliding window approach
(window\_length 10 000).

\emph{Chromosome Pathway Analysis}

In order to compare the gene composition in regions of interest, we
performed statistical gene enrichment analysis using the online version
of g:Profiler (version e102\_eg49\_p15\_7a9b4d6) with the g:SCS multiple
testing correction method applying a significance threshold of 0.05
(Raudvere et al., 2019). This tool includes a database of Gene Ontology
(GO) (Ashburner et al., 2000; Gene Ontology Consortium, 2021) and Kyoto
Encyclopedia of Genes and Genomes (KEGG) (Kanehisa \& Goto, 2000)
functional annotations based on the draft male mountain pine beetle
genome (DendPond\_male\_1.0). We used BLASTn with the gene models
identified in the final female assembly to identify the corresponding
gene IDs in this draft genome for analysis.

\emph{Comparison of linkage groups to SNP cohorts identified in previous
work}

Using the draft female genome assembly, Trevoy et al. (2019) conducted
principal component analyses and found that plateaus of high-loading
SNPs in linkage disequilibrium (LD) from a number of scaffolds were
driving clustering patterns on the first four principal component (PC)
axes; plateaus in PCs 1 and 3 were primarily related to geography, PC 2
was sex-linked, and PC 4 was much smaller and not clearly attributed to
geography or sex. To determine physical linkage and the chromosomal
locations of these SNPs, we assessed the correspondence between the
scaffolds from the draft female genome containing SNPs with the highest
loadings up to and including the plateaus in each PC shown in Trevoy et
al. (2019) and the final female assembly using BLAST+ v2.10.0 (Camacho
et al., 2009). For PCs 1 and 2, we included draft scaffolds containing
SNPs that had loadings equal to or greater than 0.05; for PC 3 this
cut-off was 0.087, and for PC 4, it was 0.1. We created a custom BLAST
database out of the final female assembly, and then used BLASTn to query
the draft scaffolds for each PC against the new assembly, specifying a
minimum e-value of 10\textsuperscript{-5}. For each PC, hits were sorted
first based on e-value and then bitscore, outputting the single best
match to the final assembly for each draft assembly scaffold.

\emph{Assessment of neo-sex chromosome architecture}

Identification of the neo-Y and ancestral-X portions of the neo-sex
chromosomes may clarify the process of sex chromosome evolution in
mountain pine beetle. Previous work has identified likely neo-Y
scaffolds (Bracewell et al., 2017; Dowle et al., 2017) and ancestral-X
scaffolds (Keeling et al., 2013c) in the draft male genome assembly. We
assessed the genomic locations of the putative neo-Y scaffolds in the
final male genome assembly with BLASTn as described above, using the
final male assembly to create a custom database and then blasting the
draft scaffolds identified by Dowle et al. (2017) against it. We also
used NUCmer to identify the regions of the neo-X chromosome derived from
the ancestral-X scaffolds of the draft male assembly predicted by
Keeling et al. (2013c).

\textbf{Results}

\emph{Genome assembly and annotation}

The original draft genomes used paired-end and mate-pair Illumina
library sequencing (Keeling et al., 2013c). We made substantial
improvements to these assemblies with proximity ligation-based
scaffolding with HiRise; linkage-map/ALLMAPS-informed corrections and
scaffolding; further improvements with LINKS, RAILS, and ABySS-Sealer;
and PhylOligo-based removal of contaminant scaffolds (Fig. 1, Tab. 1).
All of these tools were developed after the original draft assemblies
were prepared. A comparison of scaffold sizes between draft and final
genome assemblies is shown in Supp. Fig. 2. The final female and male
genome assembly sizes were 223.7 and 224.8 Mb, with N50s/L50s of 16.6
Mb/4 and 16.4 Mb/4, respectively. Gregory et al. (2013) used flow
cytometry to estimate a 208 Mb genome size. The non-N portions of the
genome assemblies were very similar to this value, 214.0 Mb for the
female assembly and 210.5 Mb for the male assembly. Compared to the
draft assemblies, N50 values increased by 26- and 36-fold, and the
number of scaffolds decreased by 67 and 75 percent, respectively. Ninety
percent of each assembly was contained in the largest 12 (female) and 11
(male) scaffolds. Based on linkage mapping information, these 12 largest
scaffolds in the female assembly represent the karyotype of this species
(11 AA + neo-XX). The male assembly did not contain a large scaffold
representing the neo-Y chromosome.

Each step in the assembly process contributed to the improved
assemblies, and incremental assembly statistics at each step are shown
in Supp. Tab. 1. Chicago HiRise scaffolding dramatically increased
contiguity, reducing the number of scaffolds by 56-66\%. Hi-C HiRise
scaffolding reduced the number of scaffolds by an additional 21\%. The
linkage map information allowed us to correct misjoins in the HiRise
assemblies and join additional scaffolds. Visualization of the linkage
map information with ALLMAPS allowed us to identify several instances
where scaffolds from the Chicago HiRise step were flipped and/or
out-of-order with adjacent scaffolds compared to the linkage map
information and the assembly from the other sex when they were
scaffolded at the Hi-C HiRise step, even though both assemblies were
based upon the same scaffolding information. An example is shown in Fig.
2. In total, nine of the twelve largest scaffolds were modified (Supp.
Fig. 3).

In one case only, a scaffold from the draft male assembly was flipped
and misplaced during the earlier Chicago HiRise step. Based on linkage
map information, ALLMAPS joined three scaffolds to make the neo-X in the
female assembly, and four scaffolds to make the neo-X in the male
assembly. This made the neo-X scaffold the largest scaffold in both
final assemblies. The LINKS scaffolding step made only two and six
joins, the RAILS step made eight and eleven joins while also filling in
18\% and 9\% of the existing gaps within scaffolds, and ABySS-Sealer
filled in 38\% and 47\% of the remaining gaps of the female and male
genomes, respectively. We then identified and removed contaminant
scaffolds with PhylOligo. These contaminant scaffolds from the female
and male assemblies matched most similarly to \emph{Serratia} spp.
and\emph{Acinetobacter} spp., respectively. Both of these genera in the
Gammaproteobacteria have been found in the bark beetle gut bacteriome
(Hern\selectlanguage{ngerman}ández-García et al., 2017). The final assemblies showed good
consistency between sexes in both shared synteny and chromosomal
arrangement (Supp. Fig. 4), and also contained 95\% of the 1367 Insecta
orthologous gene set (Insecta\_odb10, Creation date: 2020-09-10, Supp.
Fig. 5).

To annotate the genome, we used evidence from coleopteran proteins
and\emph{Dendroctonus} spp. transcripts, with \emph{ab initio} methods
for gene prediction with three rounds of Maker3. We identified 13\selectlanguage{english} 393
and 13\selectlanguage{english} 601 gene models in the female and male genomes, respectively.
This represents approximately a 4\% increase from the original draft
genome annotations. These gene models contained 91\% of the Insecta
orthologous gene set (Supp. Fig. 5) and 74\% shared significant homology
to proteins in the UniProtKB/Swiss-Prot 2020\_01 database. Repetitive
elements occupied approximately 23\% and 20\% of the female and male
genome assemblies, respectively (Supp. Tab. 2).

\emph{Linkage Mapping}

\emph{Controlled crosses}

From the wild-collected larvae, 74 females and 69 males (F0 generation)
were used to establish 66 North x South crosses, producing the F1
generation (Supp. Fig. 1). A total of 92 female and 69 male F1 beetles
were used to establish 66 full-sibling crosses, producing the F2
generation (Supp. Fig. 1). Of these full-sibling crosses, 22 (33\%)
failed to produce offspring and 10 (15\%) produced three or fewer
offspring. A further 20 (30\%) of these F1 families were discarded due
to insufficient brood size (\textless{}10 individuals) or evidence of
gallery mixing (i.e., clear family units could not be identified). Thus,
14 (21\%) F1 broods remained for producing an F2 generation. A total of
17 F1 and 212 F2 individuals were genotyped (N=229), while ten F1
specimens (nine males and one female) that were parents of the F2 broods
were not recovered from their bolts (Supp. Tab. 3).

\emph{Genotyping and Linkage Mapping}

Genotyping 229 mountain pine beetles produced 414 million reads and
22\selectlanguage{english} 242 SNPs after alignment of samples to the final female genome and
22\selectlanguage{english} 087 SNPs after alignment to the final male genome. Thirteen
individuals failed IBD testing in LepMap3, including one entire family
and three singletons from different families (Supp. Tab. 3). After
removing these individuals and imputing genotype information for the ten
missing F1 specimens, a total of 115 males and 111 females (N=226) in 13
families were used for linkage map construction.

A LOD score of 12 recovered all 11 autosomes and the neo-X chromosome in
each linkage map (Supp. Fig. 3 and Fig. 3). Increasing the LOD score to
values much higher than 12 in the male-aligned linkage maps and/or
adjusting the minimum number of SNPs required to form linkage groups
failed to recover a putative neo-Y chromosome, but did result in the
splitting of other autosomal linkage groups formed from single genomic
scaffolds. We were therefore unable to recover a neo-Y linkage group in
the male-aligned linkage map. Final female and male linkage maps
constructed from the final genome assemblies contained 2\selectlanguage{english} 795 and 2\selectlanguage{english} 910
SNPs, respectively. Linkage mapping results are shown in Supp. Tables
4-6.

The recovered linkage groups in the female-aligned linkage map were
consistent with those in the male map except for female chromosome 10,
which was equivalent to male chromosome 11 and vice versa (Fig. 3).
Genomic scaffold numbering in this study reflects assembled chromosome
size (with chromosome 1 being the largest scaffold) and is therefore not
necessarily consistent between the male and female assemblies. In this
instance, although male chromosome 11 and female chromosome 10 were
syntenous, the assembled size in base pairs of the male scaffold was
smaller than the female scaffold. The female-aligned linkage groups
ranged from approximately 50 to 90 cM in size, and the equivalent
male-aligned linkage map ranged from 50 to 70 cM. This size difference
can be attributed to the male and female neo-X linkage groups; the
female neo-X was approximately 20 cM larger than the male copy.
Chromosome 4 in both the male- and female-aligned linkage maps contained
a region approximately 30 cM in size with strong genetic linkage
(approximately 0.1 cM/locus), although the exact size and chromosomal
location of this region differed slightly between the maps (Fig. 3).

\emph{Assessment of local adaptation in Canadian mountain pine beetle
populations}

\emph{Population Differentiation and Tajima's D}

We estimated genome-wide patterns of diversity (\selectlanguage{greek}π\selectlanguage{english}), differentiation
(\emph{F\textsubscript{ST}} ) and Tajima's \emph{D} of the pool-seq data
using non-overlapping 10 kb sliding windows for northern and southern
beetle populations. Average pairwise \emph{F\textsubscript{ST}} across
the female genome was low at \emph{F\textsubscript{ST}} = 0.063 (SD =
0.038), with averages for individual chromosomes ranging from 0.053 (SD
= 0.042) to 0.092 (SD = 0.062) (Fig. 4a, Tab. 2). Average
Tajima's\emph{D} was negative for both northern (-0.859) and southern
(-0.808) populations (Fig. 4b). Two regions exhibited elevated genetic
differentiation. The first region was located at the terminal end of the
neo-X chromosome between 47 and 64 Mb (Fig.
4c).\emph{F\textsubscript{ST}} for this region was significantly higher
than the rest of the chromosome (average \emph{F\textsubscript{ST}} =
0.096 versus 0.038, unpaired \emph{t} -test: \emph{t} = 42.497, \emph{p}
-value = 5 x10\textsuperscript{-276}). Average Tajima's \emph{D} in this
genomic region was also significantly lower in the northern population
compared to the rest of the chromosome (-1.294 versus -0.913,
unpaired\emph{t} -test: \emph{t} = 20.771, \emph{p} -value = 5.3
x10\textsuperscript{-88}). This difference in Tajima's \emph{D} was not
present in the southern population.

The second region of elevated genetic differentiation was on chromosome
4 between 2.7 and 6.7 Mb (average \emph{F\textsubscript{ST}} = 0.193
versus 0.061, unpaired \emph{t} -test: \emph{t} = 64.535, \emph{p}
-value = 5.1 x10\textsuperscript{-239}) (Fig. 4d). This second region
exhibited a significant difference in Tajima's \emph{D} between northern
and southern populations (Tajima's \emph{D} = -0.615 versus 0.109,
unpaired\emph{t} -test: \emph{t} = 27.756, \emph{p} -value = 7.5
x10\textsuperscript{-118}) and was located within the region of strong
genetic linkage in female-aligned linkage group 4 (Fig. 3a). In the
linkage map, this region was roughly 23 to 30 cM in size and contained
147 of the 515 SNPs (approximately 29\%) in this linkage group. Upstream
of this second region (start of chromosome to 2.7 Mb) we detected
significantly lower levels of \emph{F\textsubscript{ST}} compared to the
rest of the chromosome located downstream of this region (6.7 Mb to end
of chromosome) (average \emph{F\textsubscript{ST}} = 0.045 versus 0.065,
unpaired \emph{t} -test: \emph{t} = 19.043, \emph{p} -value =
2.1x10\textsuperscript{-68}). The Tajima's \emph{D} in this upstream
region also showed a significant difference between northern and
southern populations (Tajima's \emph{D} = 0.076 versus -0.155,
unpaired\emph{t} -test: \emph{t} = 5.841, \emph{p} -value = 8.76
x10\textsuperscript{-9}). No differences between populations in
Tajima's\emph{D} were apparent in the downstream region.

We characterized the gene composition of these two chromosomes using
Gene Ontology (GO) annotation and KEGG pathway enrichment analysis.
Overall, the neo-X chromosome contained genes significantly enriched in
the Cellular Component GO terms: cell periphery (GO:0071944), plasma
membrane (GO:0005886), and SWI/SNF complex (GO:0016514) (Supp. File 1).
The 18 Mb region at the terminal end of the neo-X chromosome, containing
880 genes, was significantly enriched in genes with the Molecular
Function GO term: catalytic activity, acting on a protein (GO:0140096);
Biological Process GO terms: transport (GO:0006810) and establishment of
localization (GO:0051234); and the Cellular Component GO terms: SWI/SNF
complex (GO:0016514) and gap junction (GO:0005921). No KEGG pathway was
significantly enriched in this region, or the neo-X chromosome overall.

On chromosome 4, neither the region of elevated differentiation between
2.7 and 6.7 Mb, nor the region downstream had any significantly enriched
GO terms. However, the region upstream contained genes significantly
enriched in the Molecular Function GO terms: alcohol-forming fatty
acyl-CoA reductase activity (GO:0102965), fatty-acyl-CoA reductase
(alcohol-forming) activity (GO:0080019), oxidoreductase activity, acting
on the aldehyde or oxo group of donors, NAD or NADP as acceptor
(GO:0016620), and oxidoreductase activity, acting on the aldehyde or oxo
group of donors (GO:0016903). No KEGG pathways were significantly
enriched for chromosome 4.

\emph{Comparison of linkage groups to SNP cohorts identified in previous
work}

Using the PC loading plateaus described above and in Trevoy et al.
(2019), we identified and extracted scaffolds from the draft female
assembly (Keeling et al., 2013c) that had the highest loadings for each
of the first four axes (Supp. Tab. 7). The 48 draft scaffolds identified
from PC 1, which Trevoy et al. (2019) found to be related to geography,
contained 101 SNPs and had top BLAST hits to 16 scaffolds in the final
female assembly. Chromosome 4 contained 66 of these SNPs. PC 2 contained
214 SNPs on 60 draft scaffolds. With only one exception, all of these
draft scaffolds had their top BLAST hits on the neo-X chromosome in the
female assembly. This is consistent with the sex-based clustering along
PC 2 found in Trevoy et al. (2019). PC 3 contained 50 SNPs on 11 draft
scaffolds. These draft scaffolds aligned to chromosome 4, and 9 of them
were also found in PC 1. Of the 68 unique SNPs from PCs 1 and 3 that
corresponded to female chromosome 4, 13 were located within the region
of elevated \emph{F\textsubscript{ST}} and Tajima's \emph{D} on that
chromosome, as described above (Supp. Tab. 8). Finally, two draft
scaffolds containing 25 SNPs from PC 4 matched to chromosome 9 in the
final female assembly.

\emph{Assessment of neo-sex chromosome architecture}

We previously identified six large scaffolds in the draft male assembly
with low nucleotide variation that we predicted constituted the
ancestral-X portion of the neo-X chromosome (Keeling et al., 2013c). In
the new assemblies, these scaffolds were found at the terminal end of
the neo-X scaffold in the same orientation and almost the same order as
previously predicted, starting at approximately 58.3 Mb and 51.1 Mb for
the male and female neo-X scaffolds, respectively (Supp. Fig. 6).

We also used BLASTn to determine where 800 putative neo-Y scaffolds from
the draft assembly found by Dowle et al. (2017) were located in the new
male assembly. Of these draft scaffolds, 454 matched to the male neo-X
chromosome, and the remaining 346 draft scaffolds corresponded to 202
other scaffolds in the final genome (Supp. Tab. 9). The male neo-X
chromosome spans 70 670 744 bp in length, and the 454 draft scaffolds
corresponding to the neo-X chromosome were located between base pairs
1 440 820 and 59 293 732 (Supp. Tab. 10). These draft scaffolds were not
contiguous but were instead interspersed along this region of the neo-X
chromosome (Fig. 5). The remaining portion of the final neo-X chromosome
downstream of base pair 59 293 732 lacked any putative neo-Y scaffolds
and also broadly overlapped with the region of
elevated\emph{F\textsubscript{ST}} and Tajima's \emph{D} on the final
female neo-X chromosome, described above. All these data consistently
support the same location of the ancestral-X at the terminal end of the
neo-X chromosome.

\textbf{Discussion}

\emph{Genome Assembly}

The draft mountain pine beetle genomes were assembled a decade ago from
paired-end and mate-pair Illumina sequences during the early development
of long-read sequencing and proximity ligation sequencing technologies.
Since then, analytical tools have improved and new assembly tools have
been developed to incorporate long-read and proximity ligation
sequencing technologies (Amarasinghe et al., 2020). An upsurge of
chromosome-level assemblies is currently being published using these new
approaches and tools, including for many insects. We used proximity
ligation sequencing and the improved analytical tools, along with
linkage map information to enhance the draft mountain pine beetle genome
assemblies. Most of the assembly content was present in chromosome-sized
scaffolds. Linkage map information provided critical complementary
information that highlighted inconsistencies in the proximity
ligation-based assemblies. Although we saw no evidence of draft
scaffolds being incorrectly associated with each other, we found several
instances of incorrect local ordering and orientation. This typifies the
recurring challenges of this approach. However, bioinformatic tools to
reduce these errors continue to develop and improve (e.g., Nakabayashi
\& Morishita, 2020).

At present, genomes are published for only two other bark beetle species
(Curculionidae: Scolytinae): the coffee borer beetle, \emph{Hypothenemus
hampei} (Ferrari) (Vega et al., 2015), and the European spruce bark
beetle, \emph{Ips typographus} L. (Powell et al., 2020). Although the
latter was completed with long-read PacBio sequences, our use of the
complementary approaches of proximity ligation sequencing and linkage
map data for scaffolding existing draft genomes resulted in more
contiguous assemblies. However, gene content was similarly complete
between assemblies (95.6\% for \emph{I. typographus} versus 95.2\% and
94.8\%, for the female and male mountain pine beetle assemblies,
respectively, based upon BUSCO v4.1.4 with insecta\_odb10, created
2020-09-10).

\emph{Genomic architecture reflects geographic divergence and
adaptation}

Divergence between northern and southern Canadian populations of
mountain pine beetle has been shown in several population genetic
studies (Batista et al., 2016; Cullingham et al., 2011; Janes et al.,
2014). Using available genome-wide pool-seq data (Keeling et al., 2013c)
and our newly assembled genomes, we evaluated this genetic
differentiation at the chromosome level. The mountain pine beetle genome
displayed two notable regions of elevated genetic divergence. The first
divergent region is on the terminal end of the neo-X chromosome (Fig.
4c; Fig. 5a). This 18 Mb region represents the ancestral-X portion of
the neo-X chromosome. Genes in this region were significantly enriched
in GO terms associated with transport, chromatin remodeling, gap
junctions, and catalytic activity (Supp. file 1).

Previous work by Bracewell et al. (2017) on mountain pine beetle
populations in the northern US contrasts with our findings on genetic
divergence in Canadian populations, as they did not identify differences
in genetic divergence between the ancestral-autosomal and ancestral-X
regions of the neo-X chromosome. However, that study reported higher
genetic diversity in the ancestral-autosomal region. This differs from
our findings where the ancestral-X region had the lowest levels of
nucleotide diversity (Tab. 2). We also found the ancestral-autosomal
region of neo-X had relatively similar levels of diversity to those of
the other autosomes, but also showed the lowest levels of
differentiation. In each chromosome and sub-region, we saw a higher
level of genetic diversity compared to Bracewell et al. (2017), which
may be explained by the recent population expansion in Canadian
populations.

The second region with substantially increased divergence involves 4 Mb
near the beginning of chromosome 4 (Fig. 4d). This region shows a
north-south difference in Tajima's \emph{D} , with the southern
population having positive Tajima's \emph{D} values compared to the
negative-shifted values in the northern population that are typical of
most of the rest of chromosome 4 and other chromosomes. Positive
Tajima's \emph{D} values in the southern population suggest that genes
within this region are undergoing balancing selection, while negative
Tajima's \emph{D} values in the north highlight an excess of rare
alleles in this region. It remains unknown if these polymorphisms play a
role in population expansion or local adaptation of mountain pine
beetles in northern British Columbia and Alberta. It is possible that
this pattern of differentiation indicates a chromosomal inversion
between the northern and southern populations, with both inversion
variants being present in the southern population. This could alter the
selection occurring on genes within this genomic region and could also
alter patterns of gene expression near the inversion breakpoints (Durmaz
et al., 2021). However, our linkage mapping results did not indicate
patterns of recombination on chromosome 4 that were consistent with an
inversion, and so the differentiated \emph{F\textsubscript{ST}} and
Tajima's \emph{D} in this region may instead reflect the recent
expansion and establishment of mountain pine beetle populations north
and east of their historical Canadian range. Given the size of this
differentiated region and the biological importance of the genes
involved in protein synthesis and gene regulation found to be enriched
in this part of the chromosome, it presents an interesting target for
investigating differences in local adaptation between populations of the
mountain pine beetle.

\emph{Neo-sex chromosome characterization and its applications to
studies of chromosomal evolution}

Alignment of the draft neo-Y scaffolds to the neo-X revealed the
approximate boundary between the ancestral-autosomal and ancestral-X
portions of the neo-X chromosome. Our data suggest that the ancestral-X
region may be slightly larger than 10 Mb (Fig. 5), or approximately 14\%
of the total length of the neo-X chromosome.

Using linkage mapping and genome assemblies, we found that markers in
each of the first four principal components (PCs) identified in Trevoy
et al. (2019) were physically linked to three chromosomes. These
findings validate the use of PC loadings exhibiting marker plateaus to
infer genomic linkage in the absence of a linkage map or
chromosome-scale genome assembly. Trevoy et al. (2019) also showed that
high PC 2 loadings distinguished SNPs in beetles that were almost
entirely heterozygous in males and homozygous in females. They
hypothesized that this reflects sex-specific nucleotide changes in
paralogous neo-X linked versus neo-Y linked genes. If true, the
male-specific heterozygosity found by Trevoy et al. (2019) actually
resulted from erroneous alignment of neo-Y reads to the
ancestral-autosomal region of the neo-X. Support for this conclusion
comes from neo-Y scaffold locations on the draft genome, as identified
by Dowle et al. (2017), mapping to our neo-X chromosome (Fig. 5).
Further, if any of the draft genome scaffolds containing high-loading
SNPs on the PC 2 axis identified by Trevoy et al. (2019) truly
represented the male copy of the neo-X, we would expect some of them to
align to the ancestral-X region as well, but this was not the case.
Rather, we identify a mosaic-like pattern across neo-Y scaffolds that is
consistent with both the recent origin of mountain pine beetle neo-sex
chromosomes and the gradual accumulation of chromosome-specific
mutations.

Karyotype varies across \emph{Dendroctonus} and includes both typical
Xy\textsubscript{p} sex chromosomes (as found for the majority of
beetles) and neo-sex chromosomes (Z\selectlanguage{ngerman}úñiga et al., 2002), like in mountain
pine beetle. We have generated the first chromosome-level assembly in
this genus, and others are forthcoming (Casola et al., 2020; Keeling et
al., 2020). The mountain pine beetle neo-sex chromosomes represent
approximately 30\% of its genome content. Genome assemblies from the
other \emph{Dendroctonus} spp. should identify which ancestral autosomes
became part of the neo-sex chromosomes in mountain pine beetle, and the
consequences of such large sex chromosomes for mountain pine beetle
physiology and ecology.

This improved mountain pine beetle genome may provide another window for
studying chromosome evolution relating to neo-XY development and Y
chromosome degeneration. Currently, much of our understanding of neo-XY
systems comes from studies of \emph{Drosophila} spp. (Wei \& Bachtrog,
2019), while mammalian and \emph{Drosophila} systems provide insight
into Y degeneration (e.g., Muller's ratchet) (Charlesworth \&
Charlesworth, 2000). However, neither mammals nor \emph{Drosophila} are
ideal systems for these studies because divergence times among species
can be considerable -- their Y chromosomes tend to be quite old,
exhibiting substantial degeneration (Charlesworth \& Charlesworth,
2000). Thus, evolutionarily more recent neo-Y systems, like that of
mountain pine beetle, are particularly valuable.

Several other features of the mountain pine beetle genome, coupled with
the improved genomic resource presented here, make it an interesting
chromosomal evolution study system. First, compared to traditional
chiasmatic or achiasmatic XY systems, the ancestral
Xy\textsubscript{p}karyotypes of many \emph{Dendroctonus} spp. may more
readily result in neo-XY systems. (Blackmon \& Demuth, 2014; Dutrillaux
\& Dutrillaux, 2017). Second, Y turnover, due to degeneration followed
by neo-Y gains, is common in many beetle groups, and is approximately
34\% in\emph{Dendroctonus} (Blackmon \& Demuth, 2014). Third, neo-Y
haplotype groups corresponding to geography and low hybrid viability
have been identified in the mountain pine beetle (Dowle et al., 2017).
Thus, further study of the mountain pine beetle genome may allow us to
answer a variety of questions relating to Muller's ratchet-related
concepts, and implications for neo-Y chromosome turnover, in
phylogeographic and cyclic population irruption contexts.

In conclusion, we completed proximity ligation-based scaffolding of the
draft genome of \emph{D. ponderosae} , supported by linkage mapping, to
generate chromosome-level assemblies. These assemblies support
genome-level investigations of many biological processes for this
keystone species and silvicultural pest and will serve as a valuable
resource for functional and evolutionary studies of
other\emph{Dendroctonus} species, other Scolytinae and Coleoptera more
generally.

\textbf{Acknowledgements}

This research was supported by a Genomics Research and Development
Initiative (GRDI) grant from the Government of Canada (to CIK), a
Natural Sciences and Engineering Research Council of Canada (NSERC)
Discovery Grant (RGPIN-2018-04920, to FAHS), and an NSERC Strategic
Partnership Grants for Networks (NET GP 434810-12, to DPWH and FAHS) to
the TRIA Network, with contributions from Alberta Agriculture and
Forestry, fRI Research, Manitoba Conservation and Water Stewardship,
Natural Resources Canada---Canadian Forest Service, Northwest
Territories Environment and Natural Resources, Ontario Ministry of
Natural Resources and Forestry, Saskatchewan Ministry of Environment,
West Fraser Timber Co. Ltd. and Weyerhaeuser Canada Ltd. This research
was enabled in part by computing support provided by WestGrid
(www.westgrid.ca) and Compute Canada Calcul Canada
(www.computecanada.ca).

\textbf{References}

Amarasinghe, S. L., Su, S., Dong, X., Zappia, L., Ritchie, M. E., \&
Gouil, Q. (2020). Opportunities and challenges in long-read sequencing
data analysis. \emph{Genome Biology, 21} (1), 30.
https://doi.org/10.1186/s13059-020-1935-5

Amman, G. D. (1972). Mountain pine beetle brood production in relation
to thickness of lodgepole pine phloem. \emph{Journal of Economic
Entomology, 65} , 138-140.

Andersson, M. N., Grosse-Wilde, E., Keeling, C. I., Bengtsson, J. M.,
Yuen, M. M. S., Li, M., Hillbur, Y., Bohlmann, J., Hansson, B. S., \&
Schlyter, F. (2013). Antennal transcriptome analysis of the chemosensory
gene families in the tree killing bark beetles, \emph{Ips
typographus}and \emph{Dendroctonus ponderosae} (Coleoptera:
Curculionidae: Scolytinae). \emph{BMC Genomics, 14} , 198.

Andersson, M. N., Keeling, C. I., \& Mitchell, R. F. (2019). Genomic
content of chemosensory genes correlates with host range in wood-boring
beetles (\emph{Dendroctonus ponderosae} , \emph{Agrilus planipennis} ,
and\emph{Anoplophora glabripennis} ). \emph{BMC Genomics, 20} (1), 690.
https://doi.org/10.1186/s12864-019-6054-x

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H.,
Cherry, J. M., Davis, A. P., Dolinski, K., Dwight, S. S., Eppig, J. T.,
Harris, M. A., Hill, D. P., Issel-Tarver, L., Kasarskis, A., Lewis, S.,
Matese, J. C., Richardson, J. E., Ringwald, M., Rubin, G. M., \&
Sherlock, G. (2000). Gene ontology: tool for the unification of biology.
The Gene Ontology Consortium. \emph{Nature Genetics, 25} (1), 25-29.
https://doi.org/10.1038/75556

Aw, T., Schlauch, K., Keeling, C. I., Young, S., Bearfield, J. C.,
Blomquist, G. J., \& Tittiger, C. (2010). Functional genomics of
mountain pine beetle (\emph{Dendroctonus ponderosae} ) midguts and fat
bodies. \emph{BMC Genomics, 11} (1), 215.

Barth, J. M. I., Villegas\selectlanguage{english}-Rios, D., Freitas, C., Moland, E., Star, B.,
Andre, C., Knutsen, H., Bradbury, I., Dierking, J., Petereit, C.,
Righton, D., Metcalfe, J., Jakobsen, K. S., Olsen, E. M., \& Jentoft, S.
(2019). Disentangling structural genomic and behavioural barriers in a
sea of connectivity. \emph{Molecular Ecology, 28} (6), 1394-1411.
https://doi.org/10.1111/mec.15010

Bartholome, J., Mandrou, E., Mabiala, A., Jenkins, J., Nabihoudine, I.,
Klopp, C., Schmutz, J., Plomion, C., \& Gion, J.-M. (2014).
High-resolution genetic maps of \emph{Eucalyptus}
improve\emph{Eucalyptus grandis} genome assembly. \emph{New Phytologist,
206} (4), 1283-1296. https://doi.org/10.1111/nph.13150

Batista, P. D., Janes, J. K., Boone, C. K., Murray, B. W., \& Sperling,
F. A. H. (2016). Adaptive and neutral markers both show continent-wide
population structure of mountain pine beetle (\emph{Dendroctonus
ponderosae} ). \emph{Ecology and Evolution, 6} (17), 6292-6300.
https://doi.org/10.1002/ece3.2367

Bay, R. A., \& Ruegg, K. (2017). Genomic islands of divergence or
opportunities for introgression? \emph{Proceedings of the Royal Society
B: Biological Sciences, 284} (1850), 20162414.
https://doi.org/10.1098/rspb.2016.2414

Berlin, S., Hallingbck, H. R., Beyer, F., Nordh, N.-E., Weih, M., \&
Rnnberg-Wstljung, A.-C. (2017). Genetics of phenotypic plasticity and
biomass traits in hybrid willows across contrasting environments and
years. \emph{Annals of Botany, 120} (1), 87-100.
https://doi.org/10.1093/aob/mcx029

Blackmon, H., \& Demuth, J. P. (2014). Estimating tempo and mode of Y
chromosome turnover: explaining Y chromosome loss with the fragile Y
hypothesis. \emph{Genetics, 197} (2), 561-572.
https://doi.org/10.1534/genetics.114.164269

Blackmon, H., Ross, L., \& Bachtrog, D. (2016). Sex determination, sex
chromosomes, and karyotype evolution in insects. \emph{Journal of
Heredity, 108} (1), 78-93. https://doi.org/10.1093/jhered/esw047

Bleidorn, C. (2015). Third generation sequencing: technology and its
potential impact on evolutionary biodiversity research.\emph{Systematics
and Biodiversity, 14} (1), 1-8.
https://doi.org/10.1080/14772000.2015.1099575

Bracewell, R. R., Bentz, B. J., Sullivan, B. T., \& Good, J. M. (2017).
Rapid neo-sex chromosome evolution and incipient speciation in a major
forest pest. \emph{Nature Communications, 8} , 1-14.
https://doi.org/10.1038/s41467-017-01761-4

Bushnell, B. (2020). \emph{BBMap short read
aligner}https://sourceforge.net/projects/bbmap/

Butler, J. B., Vaillancourt, R. E., Potts, B. M., Lee, D. J., King, G.
J., Baten, A., Shepherd, M., \& Freeman, J. S. (2017). Comparative
genomics of \emph{Eucalyptus} and \emph{Corymbia} reveals low rates of
genome structural rearrangement. \emph{BMC Genomics, 18} , 397.
https://doi.org/10.1186/s12864-017-3782-7

Cairns, J., Freire-Pritchett, P., Wingett, S. W., Vrnai, C., Dimond, A.,
Plagnol, V., Zerbino, D., Schoenfelder, S., Javierre, B.-M., Osborne,
C., Fraser, P., \& Spivakov, M. (2016). CHiCAGO: robust detection of DNA
looping interactions in Capture Hi-C data. \emph{Genome Biology, 17}
(1), 1-17. https://doi.org/10.1186/s13059-016-0992-2

Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J.,
Bealer, K., \& Madden, T. L. (2009). BLAST+: architecture and
applications. \emph{BMC Bioinformatics, 10} , 421.
https://doi.org/10.1186/1471-2105-10-421

Casola, C., Landa, S., Hjelman, C., Jonika, M., Sullivan, B., Kyre, B.,
Rieske-Kinney, L. K., \& Blackmon, H. (2020). \emph{Using comparative
genomics to identify species-specific targets for RNAi in Dendroctonus
bark beetles} 2020 Entomology Virtual Annual Meeting of the
Entomological Society of America,

Catchen, J. M., Amores, A., Hohenlohe, P., Cresko, W., \& Postlethwait,
J. H. (2011). Stacks: building and genotyping loci \emph{de novo} from
short-read sequences. \emph{G3: Genes Genomes Genetics, 1} , 171-182.
https://doi.org/10.1534/g3.111.000240/-/DC1

Charlesworth, B., \& Charlesworth, D. (2000). The degeneration of Y
chromosomes. \emph{Philosophical Transactions Of The Royal Society Of
London B Biological Sciences, 355} .

Chiu, C. C., Keeling, C. I., \& Bohlmann, J. (2018). Monoterpenyl esters
in juvenile mountain pine beetle and sex-specific release of the
aggregation pheromone \emph{trans} -verbenol. \emph{Proceedings of the
National Academy of Science, USA, 115} (14), 3652-3657.
https://doi.org/10.1073/pnas.1722380115

Chiu, C. C., Keeling, C. I., \& Bohlmann, J. (2019a). The cytochrome
P450 CYP6DE1 catalyzes the conversion of \selectlanguage{greek}α-\selectlanguage{english}pinene into the mountain pine
beetle aggregation pheromone \emph{trans} -verbenol. \emph{Scientific
Reports, 9} (1), 1477.

Chiu, C. C., Keeling, C. I., \& Bohlmann, J. (2019b). Cytochromes P450
preferentially expressed in antennae of the mountain pine
beetle.\emph{Journal of Chemical Ecology, 45} (2), 178-186.
https://doi.org/10.1007/s10886-018-0999-0

Chiu, C. C., Keeling, C. I., \& Bohlmann, J. (2019c). Functions of
mountain pine beetle cytochromes P450 CYP6DJ1, CYP6BW1 and CYP6BW3 in
the oxidation of pine monoterpenes and diterpene resin acids. \emph{PLoS
ONE, 15} (5), e0216753.

Christmas, M. J., Wallberg, A., Bunikis, I., Olsson, A., Wallerman, O.,
\& Webster, M. T. (2019). Chromosomal inversions associated with
environmental adaptation in honeybees. \emph{Molecular Ecology, 28} (6),
1358-1374. https://doi.org/10.1111/mec.14944

Chu, J. (2020). \emph{JupiterPlot} .
https://github.com/JustinChu/JupiterPlot

Cullingham, C. I., Cooke, J. E., Dang, S., Davis, C. S., Cooke, B. J.,
\& Coltman, D. W. (2011). Mountain pine beetle host-range expansion
threatens the boreal forest. \emph{Molecular Ecology, 20} (10),
2157-2171. https://doi.org/10.1111/j.1365-294X.2011.05086.x

Cullingham, C. I., Janes, J. K., Hamelin, R. C., James, P. M. A.,
Murray, B. W., \& Sperling, F. A. H. (2018). The contribution of
genetics and genomics to understanding the ecology of the mountain pine
beetle system. \emph{Canadian Journal of Forest Research, 49} (7),
721-730. https://doi.org/10.1139/cjfr-2018-0303

Dfam Consortium. (2019). \emph{Dfam TE Tools Container} . Retrieved Dec.
17, 2019 from https://github.com/Dfam-consortium/TETools

Dowle, E. J., Bracewell, R. R., Pfrender, M. E., Mock, K. E., Bentz, B.
J., \& Ragland, G. J. (2017). Reproductive isolation and environmental
adaptation shape the phylogeography of mountain pine beetle
(\emph{Dendroctonus ponderosae} ). \emph{Molecular Ecology, 26} (21),
6071-6084. https://doi.org/10.1111/mec.14342

Durmaz, E., Kerdaffrec, E., Katsianis, G., Kapun, M., \& Flatt, T.
(2021). How selection acts on chromosomal inversions. In John Wiley \&
Sons, Ltd (Ed.), \emph{eLS} . John Wiley \& Sons, Ltd,.
https://doi.org/doi.org/10.1002/9780470015902.a0028745

Dutrillaux, A. M., \& Dutrillaux, B. (2017). Evolution of the sex
chromosomes in beetles. I. The loss of the Y
chromosome.\emph{Cytogenetic and Genome Research, 152} (2), 97-104.
https://doi.org/10.1159/000478075

Ferretti, L., Ramos-Onsins, S. E., \& Perez-Enciso, M. (2013).
Population genomics from pool sequencing. \emph{Molecular Ecology, 22}
(22), 5561-5576. https://doi.org/10.1111/mec.12522

Fierst, J. L. (2015). Using linkage maps to correct and scaffold\emph{de
novo} genome assemblies: methods, challenges, and computational tools.
\emph{Frontiers in Genetics, 6} , 1-8.
https://doi.org/10.3389/fgene.2015.00220

Fonseca, J. (2019). \emph{xdot: Interactive viewer for Graphviz dot
files.} In (Version 1.1) https://pypi.org/project/xdot/

Gagnaire, P. A., Pavey, S. A., Normandeau, E., \& Bernatchez, L. (2013).
The genetic architecture of reproductive isolation during
speciation-with-gene-flow in lake whitefish species pairs assessed by
RAD sequencing. \emph{Evolution, 67} (9), 2483-2497.
https://doi.org/10.1111/evo.12075

Gene Ontology Consortium. (2021). The Gene Ontology resource: enriching
a GOld mine. \emph{Nucleic Acids Research, 49} (D1), D325-D334.
https://doi.org/10.1093/nar/gkaa1113

Goodsman, D. W., Koch, D., Whitehouse, C., Evenden, M. L., Cooke, B. J.,
\& Lewis, M. A. (2016). Aggregation and a strong Allee effect in a
cooperative outbreak insect. \emph{Ecological Applications, 26} (8),
2623-2636. https://doi.org/10.1002/eap.1404

Gregory, T. R., Nathwani, P., Bonnett, T. R., \& Huber, D. P. (2013).
Sizing up arthropod genomes: an evaluation of the impact of
environmental variation on genome size estimates by flow cytometry and
the use of qPCR as a method of estimation. \emph{Genome, 56} (9),
505-510. https://doi.org/10.1139/gen-2013-0044

Hern\selectlanguage{ngerman}ández-García, J. A., Briones-Roblero, C. I., Rivera-Orduña, F. N.,
\& Zúñiga, G. (2017). Revealing the gut bacteriome of\emph{Dendroctonus}
bark beetles (Curculionidae: Scolytinae): diversity, core members and
co-evolutionary patterns. \emph{Scientific Reports, 7} (1), 13864.
https://doi.org/10.1038/s41598-017-14031-6

Holt, C., \& Yandell, M. (2011). MAKER2: an annotation pipeline and
genome-database management tool for second-generation genome
projects.\emph{BMC Bioinformatics, 12} , 491.
https://doi.org/10.1186/1471-2105-12-491

James, P. M., Coltman, D. W., Murray, B. W., Hamelin, R. C., \&
Sperling, F. A. H. (2011). Spatial genetic structure of a symbiotic
beetle-fungal system: toward multi-taxa integrated landscape
genetics.\emph{PLoS ONE, 6} (10), e25359.
https://doi.org/10.1371/journal.pone.0025359

James, P. M. A., Janes, J. K., Roe, A. D., \& Cooke, B. J. (2016).
Modeling landscape-level spatial variation in sex ratio skew in the
mountain pine beetle (Coleoptera: Curculionidae). \emph{Environmental
Entomology, 45} (4), 790-801. https://doi.org/10.1093/ee/nvw048

Janes, J. K., Li, Y., Keeling, C. I., Yuen, M. M. S., Boone, C. K.,
Cooke, J. E. K., Bohlmann, J., Huber, D. P. W., Murray, B. W., Coltman,
D. W., \& Sperling, F. A. H. (2014). How the mountain pine beetle
(\emph{Dendroctonus ponderosae} ) breached the Canadian Rocky
Mountains.\emph{Molecular Biology and Evolution, 31} (7), 1803-1815.
https://doi.org/10.1093/molbev/msu135

Janes, J. K., Roe, A. D., Rice, A. V., Gorrell, J. C., Coltman, D. W.,
Langor, D. W., \& Sperling, F. A. H. (2016). Polygamy and an absence of
fine-scale structure in \emph{Dendroctonus ponderosae} (Hopk.)
(Coleoptera: Curcilionidae) confirmed using molecular
markers.\emph{Heredity, 116} , 68-74.
https://doi.org/10.1038/hdy.2015.71

Jones, P., Binns, D., Chang, H. Y., Fraser, M., Li, W., McAnulla, C.,
McWilliam, H., Maslen, J., Mitchell, A., Nuka, G., Pesseat, S., Quinn,
A. F., Sangrador-Vegas, A., Scheremetjew, M., Yong, S. Y., Lopez, R., \&
Hunter, S. (2014). InterProScan 5: genome-scale protein function
classification. \emph{Bioinformatics, 30} (9), 1236-1240.
https://doi.org/10.1093/bioinformatics/btu031

Kanehisa, M., \& Goto, S. (2000). KEGG: Kyoto encyclopedia of genes and
genomes. \emph{Nucleic Acids Research, 28} (1), 27-30.
https://doi.org/10.1093/nar/28.1.27

Keeling, C. I., Campbell, E. O., Willsey, T., Strong, W., Huber, D. P.
W., Sperling, F., \& Bleiker, K. (2020). \emph{Genome sequencing in bark
beetles} 2020 Entomology Virtual Annual Meeting of the Entomological
Society of America,

Keeling, C. I., Chiu, C. C., Aw, T., Li, M., Henderson, H., Tittiger,
C., Weng, H. B., Blomquist, G. J., \& Bohlmann, J. (2013a). Frontalin
pheromone biosynthesis in the mountain pine beetle, \emph{Dendroctonus
ponderosae} , and the role of isoprenyl diphosphate
synthases.\emph{Proceedings of the National Academy of Science, USA,
110} (47), 18838-18843. https://doi.org/10.1073/pnas.1316498110

Keeling, C. I., Henderson, H., Li, M., Dullat, H. K., Ohnishi, T., \&
Bohlmann, J. (2013b). CYP345E2, an antenna-specific cytochrome P450 from
the mountain pine beetle, \emph{Dendroctonus ponderosae} Hopkins,
catalyses the oxidation of pine host monoterpene volatiles. \emph{Insect
Biochemistry and Molecular Biology, 43} (12), 1142-1151.
https://doi.org/10.1016/j.ibmb.2013.10.001

Keeling, C. I., Henderson, H., Li, M., Yuen, M., Clark, E. L., Fraser,
J. D., Huber, D. P. W., Liao, N. Y., Docking, T. R., Birol, I., Chan, S.
K., Taylor, G. A., Palmquist, D., Jones, S. J. M., \& Bohlmann, J.
(2012). Transcriptome and full-length cDNA resources for the mountain
pine beetle, \emph{Dendroctonus ponderosae} Hopkins, a major insect pest
of pine forests. \emph{Insect Biochemistry and Molecular Biology, 42}
(8), 525--536.
https://doi.org/http://dx.doi.org/10.1016/j.ibmb.2012.03.010

Keeling, C. I., Li, M., Dullat, H. K., Henderson, H., Yuen, M. M. S., \&
Bohlmann, J. (2016). Quantitative metabolome, proteome and transcriptome
analysis of midgut and fat body tissues in the mountain pine
beetle,\emph{Dendroctonus ponderosae} Hopkins, and insights into
pheromone biosynthesis. \emph{Insect Biochemistry and Molecular Biology,
70} , 170-183.

Keeling, C. I., Yuen, M. M. S., Liao, N. Y., Docking, T. R., Chan, S.
K., Taylor, G. A., Palmquist, D. L., Jackman, S. D., Nguyen, A., Li, M.,
Henderson, H., Janes, J. K., Zhao, Y., Pandoh, P., Moore, R., Sperling,
F. A. H., Huber, D. P. W., Birol, I., Jones, S. J. M., \& Bohlmann, J.
(2013c). Draft genome of the mountain pine beetle, \emph{Dendroctonus
ponderosae} Hopkins, a major forest pest. \emph{Genome Biology, 14} ,
R27.

Kofler, R., Pandey, R. V., \& Schlotterer, C. (2011). PoPoolation2:
identifying differentiation between populations using sequencing of
pooled DNA samples (Pool-Seq). \emph{Bioinformatics, 27} (24),
3435-3436. https://doi.org/10.1093/bioinformatics/btr589

Korf, I. (2004). Gene finding in novel genomes {[}Research Support,
Non-U.S. Gov't

Research Support, U.S. Gov't, P.H.S.{]}. \emph{BMC Bioinformatics, 5} ,
59. https://doi.org/10.1186/1471-2105-5-59

Kriventseva, E. V., Kuznetsov, D., Tegenfeldt, F., Manni, M., Dias, R.,
Simão, F. A., \& Zdobnov, E. M. (2019). OrthoDB v10: sampling the
diversity of animal, plant, fungal, protist, bacterial and viral genomes
for evolutionary and functional annotations of orthologs. \emph{Nucleic
Acids Research, 47} (D1), D807-D811. https://doi.org/10.1093/nar/gky1053

Krzywinski, M., Schein, J., Birol, I., Connors, J., Gascoyne, R.,
Horsman, D., Jones, S. J., \& Marra, M. A. (2009). Circos: an
information aesthetic for comparative genomics {[}Research Support,
Non-U.S. Gov't{]}. \emph{Genome Research, 19} (9), 1639-1645.
https://doi.org/10.1101/gr.092759.109

Kurtz, S., Phillippy, A., Delcher, A. L., Smoot, M., Shumway, M.,
Antonescu, C., \& Salzberg, S. L. (2004). Versatile and open software
for comparing large genomes. \emph{Genome Biology, 5} (2), R12.
https://doi.org/10.1186/gb-2004-5-2-r12

Lanier, G. N. (1981). Cytotaxonomy of \emph{Dendroctonus} . In M. W.
Stock (Ed.), \emph{Application of Genetics and Cytology in Insect
Systematics and Evolution, Proceedings of the 1980 Annual Meeting of the
Entomological Society of America} (pp. 33-66). Forest, Wildlife and
Range Experimental Station, University of Idaho, Moscow.

Lanier, G. N., \& Wood, D. L. (1968). Controlled mating, karyology,
morphology, and sex-ratio in the \emph{Dendroctonus ponderosae}
complex.\emph{Annals of the Entomological Society of America, 61} (2),
517-526.

Li, C., Lin, F., An, D., Wang, W., \& Huang, R. (2018). Genome
sequencing and assembly by long reads in plants. \emph{Genes, 9} (1),
6-14. https://doi.org/10.3390/genes9010006

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N.,
Marth, G., Abecasis, G., \& Durbin, R. (2009). The Sequence
Alignment/Map format and SAMtools. \emph{Bioinformatics, 25} (16),
2078-2079. https://doi.org/10.1093/bioinformatics/btp352

Lieberman-Aiden, E., van Berkum, N. L., Williams, L., Imakaev, M.,
Ragoczy, T., Telling, A., Amit, I., Lajoie, B. R., Sabo, P. J.,
Dorschner, M. O., Sandstrom, R., Bernstein, B., Bender, M. A., Groudine,
M., Gnirke, A., Stamatoyannopoulos, J., Mirny, L. A., Lander, E. S., \&
Dekker, J. (2009). Comprehensive mapping of long-range interactions
reveals folding principles of the human genome. \emph{Science, 326}
(5950), 289-293. https://doi.org/10.1126/science.1181369

Lomsadze, A., Ter-Hovhannisyan, V., Chernoff, Y. O., \& Borodovsky, M.
(2005). Gene identification in novel eukaryotic genomes by self-training
algorithm. \emph{Nucleic Acids Research, 33} (20), 6494-6506.
https://doi.org/10.1093/nar/gki937

Ma, T., Wang, K., Hu, Q., Xi, Z., Wan, D., Wang, Q., Feng, J., Jiang,
D., Ahani, H., Abbott, R. J., Lascoux, M., Nevo, E., \& Liu, J. (2018).
Ancient polymorphisms and divergence hitchhiking contribute to genomic
islands of divergence within a poplar species complex. \emph{Proceedings
of the National Academy of Sciences of the United States of America,
115} (2), E236-E243. https://doi.org/10.1073/pnas.1713288114

Mallet, L., Bitard-Feildel, T., Cerutti, F., \& Chiapello, H. (2017).
PhylOligo: a package to identify contaminant or untargeted organism
sequences in genome assemblies. \emph{Bioinformatics, 33} (20),
3283-3285. https://doi.org/10.1093/bioinformatics/btx396

Martin, M. (2011). Cutadapt removes adapter sequences from
high-throughput sequencing reads. \emph{EMBNet Journal, 17} , 10-12.

Nadeau, J. A., Petereit, J., Tillett, R. L., Jung, K., Fotoohi, M.,
MacLean, M., Young, S., Schlauch, K., Blomquist, G. J., \& Tittiger, C.
(2017). Comparative transcriptomics of mountain pine beetle
pheromone-biosynthetic tissues and functional analysis of
CYP6DE3.\emph{BMC Genomics, 18} (1), 311.
https://doi.org/10.1186/s12864-017-3696-4

Nakabayashi, R., \& Morishita, S. (2020). HiC-Hiker: a probabilistic
model to determine contig orientation in chromosome-length scaffolds
with Hi-C. \emph{Bioinformatics, 36} (13), 3966-3974.
https://doi.org/10.1093/bioinformatics/btaa288

Ouellette, L. A., Reid, R. W., Blanchard, S. G., \& Brouwer, C. R.
(2018). LinkageMapView-rendering high-resolution linkage and QTL
maps.\emph{Bioinformatics, 34} (2), 306-307.
https://doi.org/10.1093/bioinformatics/btx576

Paulino, D., Warren, R. L., Vandervalk, B. P., Raymond, A., Jackman, S.
D., \& Birol, I. (2015). Sealer: a scalable gap-closing application for
finishing draft genomes. \emph{BMC Bioinformatics, 16} , 230.
https://doi.org/10.1186/s12859-015-0663-4

Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., \& Hoekstra,
H. E. (2012). Double digest RADseq: an inexpensive method for \emph{de
novo} SNP discovery and genotyping in model and non-model
species.\emph{PLoS ONE, 7} (5), e37135.
https://doi.org/10.1371/journal.pone.0037135

Powell, D., Gro\selectlanguage{greek}β\selectlanguage{english}e-Wilde, E., Krokene, P., Roy, A., Chakraborty, A.,
L\selectlanguage{ngerman}öfstedt, C., Vogel, H., Andersson, M. N., \& Schlyter, F. (2020). A
highly contiguous genome assembly of a major forest pest, the Eurasian
spruce bark beetle \emph{Ips typographus} . \emph{bioRxiv} .
https://doi.org/10.1101/2020.11.28.401976

Putnam, N. H., O'Connell, B. L., Stites, J. C., Rice, B. J., Blanchette,
M., Calef, R., Troll, C. J., Fields, A., Hartley, P. D., Sugnet, C. W.,
Haussler, D., Rokhsar, D. S., \& Green, R. E. (2016). Chromosome-scale
shotgun assembly using an \emph{in vitro} method for long-range
linkage.\emph{Genome Research, 26} (3), 342-350.
https://doi.org/10.1101/gr.193474.115

R Core Team. (2020). \emph{R: A language and environment for statistical
computing.} In R Foundation for Statistical Computing.
http://www.R-project.org/

Raffa, K. F., Aukema, B. H., Bentz, B. J., Carroll, A. L., Hicke, J. A.,
Turner, M. G., \& Romme, W. H. (2008). Cross-scale drivers of natural
disturbances prone to anthropogenic amplification: the dynamics of bark
beetle eruptions. \emph{Bioscience, 58} (6), 501-517.

Rastas, P. (2017). Lep-MAP3: robust linkage mapping even for
low-coverage whole genome sequencing data. \emph{Bioinformatics, 33}
(23), 3726-3732. https://doi.org/10.1093/bioinformatics/btx494

Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson,
H., \& Vilo, J. (2019). g:Profiler: a web server for functional
enrichment analysis and conversions of gene lists (2019
update).\emph{Nucleic Acids Research, 47} (W1), W191-W198.
https://doi.org/10.1093/nar/gkz369

Reid, R. W. (1962). Biology of the mountain pine
beetle,\emph{Dendroctonus monticolae} Hopkins, in east Kootenay region
of British Columbia II. Behaviour in the host, fecundity, and internal
changes in the female. \emph{Canadian Entomologist, 94} , 605-613.

Renaut, S., Grassa, C. J., Yeaman, S., Moyers, B. T., Lai, Z., Kane, N.
C., Bowers, J. E., Burke, J. M., \& Rieseberg, L. H. (2013). Genomic
islands of divergence are not affected by geography of speciation in
sunflowers. \emph{Nature Communications, 4} , 1827.
https://doi.org/10.1038/ncomms2833

Rhoads, A., \& Au, K. F. (2015). PacBio sequencing and its
applications.\emph{Genomics, Proteomics \& Bioinformatics, 13} (5),
278-289. https://doi.org/10.1016/j.gpb.2015.08.002

Robert, J. A., Bonnett, T., Pitt, C., Spooner, L. J., Fraser, J., Yuen,
M. M. S., Keeling, C. I., Bohlmann, J., \& Huber, D. P. W. (2016). Gene
expression analysis of overwintering mountain pine beetle larvae
suggests multiple systems involved in overwintering stress, cold
hardiness, and preparation for spring development. \emph{PeerJ, 4} ,
e2109. https://doi.org/10.7717/peerj.2109

Robert, J. A., Pitt, C., Bonnett, T. R., Yuen, M. M. S., Keeling, C. I.,
Bohlmann, J., \& Huber, D. P. W. (2013). Disentangling detoxification:
gene expression analysis of feeding mountain pine beetle illuminates
molecular-level host chemical defense detoxification
mechanisms.\emph{PLoS ONE, 8} (11), e77777.
https://doi.org/10.1371/journal.pone.0077777

Rochette, N. C., Rivera-Colon, A. G., \& Catchen, J. M. (2019). Stacks
2: Analytical methods for paired-end sequencing improve RADseq-based
population genomics. \emph{Molecular Ecology, 28} (21), 4737-4754.
https://doi.org/10.1111/mec.15253

Rosenberger, D. W., Venette, R. C., \& Aukema, B. H. (2016). Sexing live
mountain pine beetles \emph{Dendroctonus ponderosae} : refinement of a
behavioral method for \emph{Dendroctonus} spp. \emph{Entomologia
Experimentalis et Applicata, 160} (2), 195-199.
https://doi.org/10.1111/eea.12463

Saab, V. A., Latif, Q. S., Rowland, M. M., Johnson, T. N., Chalfoun, A.
D., Buskirk, S. W., Heyward, J. E., \& Dresser, M. A. (2014). Ecological
consequences of mountain pine beetle outbreaks for wildlife in western
North American forests. \emph{Forest Science, 60} (3), 539-559.
https://doi.org/10.5849/forsci.13-022

Safranyik, L., \& Carroll, A. L. (2006). The biology and epidemiology of
the mountain pine beetle in lodgepole pine forests. In L. Safranyik \&
B. Wilson (Eds.), \emph{The mountain pine beetle-A synthesis of biology,
management, and impacts on lodgepole pine} (pp. 3-66). Natural Resources
Canada, Canadian Forest Service.

Samarasekera, G. D. N. G., Bartell, N. V., Lindgren, B. S., Cooke, J.
E., Davis, C. S., James, P. M., Coltman, D. W., Mock, K. E., \& Murray,
B. W. (2012). Spatial genetic structure of the mountain pine beetle
(\emph{Dendroctonus ponderosae} ) outbreak in western Canada: historical
patterns and contemporary dispersal. \emph{Molecular Ecology, 21} (12),
2931-2948. https://doi.org/10.1111/j.1365-294X.2012.05587.x

Sambaraju, K. R., \& Goodsman, D. W. (2021). Mountain pine beetle: an
example of a climate-driven eruptive insect impacting conifer forest
ecosystems. \emph{CAB Reviews, 16} (018), 1-18.

Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V., \&
Zdobnov, E. M. (2015). BUSCO: assessing genome assembly and annotation
completeness with single-copy orthologs. \emph{Bioinformatics, 31} (19),
3210-3212. https://doi.org/10.1093/bioinformatics/btv351

Slater, G. S., \& Birney, E. (2005). Automated generation of heuristics
for biological sequence comparison. \emph{BMC Bioinformatics, 6} , 31.
https://doi.org/10.1186/1471-2105-6-31

Smit, A. F. A., Hubley, R., \& Green, P. (2020). \emph{RepeatMasker} .
Retrieved Sept. 15, 2018 from http://www.repeatmasker.org/

Stanke, M., Keller, O., Gunduz, I., Hayes, A., Waack, S., \&
Morgenstern, B. (2006). AUGUSTUS: \emph{ab initio} prediction of
alternative transcripts. \emph{Nucleic Acids Research, 34} , W435-W439.

Tang, H., Zhang, X., Miao, C., Zhang, J., Ming, R., Schnable, J. C.,
Schnable, P. S., Lyons, E., \& Lu, J. (2015). ALLMAPS: robust scaffold
ordering based on multiple maps. \emph{Genome Biology, 16} , 3.
https://doi.org/10.1186/s13059-014-0573-1

Tavares, H., Whibley, A., Field, D. L., Bradley, D., Couchman, M.,
Copsey, L., Elleouet, J., Burrus, M., Andalo, C., Li, M., Li, Q., Xue,
Y., Rebocho, A. B., Barton, N. H., \& Coen, E. (2018). Selection and
gene flow shape genomic islands that control floral
guides.\emph{Proceedings of the National Academy of Sciences of the
United States of America, 115} (43), 11006-11011.
https://doi.org/10.1073/pnas.1801832115

Trevoy, S. A. L., Janes, J. K., Muirhead, K., \& Sperling, F. A. H.
(2019). Repurposing population genetics data to discern genomic
architecture: A case study of linkage cohort detection in mountain pine
beetle (\emph{Dendroctonus ponderosae} ). \emph{Ecology and Evolution,
9} (3), 1147-1159. https://doi.org/10.1002/ece3.4803

van Berkum, N. L., Lieberman-Aiden, E., Williams, L., Imakaev, M.,
Gnirke, A., Mirny, L. A., Dekker, J., \& Lander, E. S. (2010). Hi-C: A
method to study the three-dimensional architecture of
genomes.\emph{Journal of Visualized Experiments} (39), 1-7.
https://doi.org/10.3791/1869

Vega, F. E., Brown, S. M., Chen, H., Shen, E., Nair, M. B.,
Ceja-Navarro, J. A., Brodie, E. L., Infante, F., Dowd, P. F., \& Pain,
A. (2015). Draft genome of the most devastating insect pest of coffee
worldwide: the coffee berry borer, \emph{Hypothenemus hampei}
.\emph{Scientific Reports, 5} , 12525. https://doi.org/10.1038/srep12525

Warren, R. L. (2016). RAILS and Cobbler: Scaffolding and automated
finishing of draft genomes using long DNA sequences. \emph{The Journal
of Open Source Software, 1} (7), 116. https://doi.org/doi:
10.21105/joss.00116

Warren, R. L., Yang, C., Vandervalk, B. P., Behsaz, B., Lagman, A.,
Jones, S. J., \& Birol, I. (2015). LINKS: Scalable, alignment-free
scaffolding of draft genomes with long reads. \emph{GigaScience, 4} ,
35. https://doi.org/10.1186/s13742-015-0076-3

Wei, K. H., \& Bachtrog, D. (2019). Ancestral male recombination
in\emph{Drosophila albomicans} produced geographically restricted neo-Y
chromosome haplotypes varying in age and onset of decay. \emph{PLoS
Genetics, 15} (11), e1008502.
https://doi.org/10.1371/journal.pgen.1008502

Zaharia, M., Bolosky, W. J., Curtis, K., Fox, A., Patterson, D.,
Shenker, S., Stoica, I., Karp, R. M., \& Sittler, T. (2011). Faster and
more accurate sequence alignment with SNAP. \emph{arXiv} , 1111.5572.
https://doi.org/arXiv:1111.5572

Zúñiga, G., Cisneros, R., Hayes, J. L., \& Macias-Samano, J. (2002).
Karyology, geographic distribution, and origin of the
genus\emph{Dendroctonus} Erichson (Coleoptera: Scolytidae). \emph{Annals
of the Entomological Society of America, 95} (3), 267-275.



\textbf{Data Accessibility}

The data that supports the findings of this study are available in the
supplementary material of this article or are openly available in NCBI
at https://www.ncbi.nlm.nih.gov with the following accession numbers.

Genome assemblies: BioProjects PRJNA638274 (male) and PRJNA638278
(female); Genome accessions JAFETF000000000 (male) and JAFETG000000000
(female).

Chicago sequencing: BioProject PRJNA638289, BioSample SAMN14918906, SRA
SRR11965410.

Hi-C sequencing: BioProject PRJNA638296, BioSample SAMN14918916, SRA
SRR11965415.

ddRAD sequencing: BioProject PRJNA694171, BioSamples
SAMN17498181-SAMN17498396.

\textbf{Author Contributions}

Designed research: CIK, EOC, PDB, VAS, SALT, DPWH, JKJ, FAHS

Performed research: CIK, EOC, PDB, VAS, SALT

Contributed new reagents or analytical tools: SALT

Analyzed data: CIK, EOC, PDB, VAS, SALT, DPWH, JKJ, FAHS

Wrote the paper: CIK, EOC, PDB, JKJ

All authors have read and approved the final version of this manuscript.



\textbf{Table 1: Summary statistics for the genome assemblies}

Summary statistics for the starting draft assemblies and the final
assemblies. Incremental summary statistics for each step in the assembly
process are shown in Supp. Tab. 1.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
\textbf{Female} & \textbf{Draft} & \textbf{Final}\tabularnewline
\midrule
\endhead
\textbf{Total size (Mb)} & 261.3 & 223.7\tabularnewline
\textbf{Non-N length (Mb)} & 212.8 & 214.0\tabularnewline
\textbf{Scaffolds:} & &\tabularnewline
\textbf{Number} & 6 520 & 2 136\tabularnewline
\textbf{Longest (Mb)} & 7.21 & 63.6\tabularnewline
\textbf{N50 (Mb)/L50} & 0.466/137 & 16.6/4\tabularnewline
\textbf{N90 (Mb)/L90} & 0.0271/943 & 3.22/12\tabularnewline
\textbf{Average length} & 40 082 & 104 747\tabularnewline
\textbf{Number \textgreater{}10 kb} & 1 696 & 554\tabularnewline
\textbf{\% of genome \textgreater{}10 kb} & 94.8 & 97.9\tabularnewline
\textbf{Gene models} & 12 873 & 13 393\tabularnewline
\textbf{Transcripts} & 13 066 & 16 041\tabularnewline
\textbf{Male} & \textbf{Draft} & \textbf{Final}\tabularnewline
\textbf{Total size (Mb)} & 252.8 & 224.8\tabularnewline
\textbf{Non-N length (Mb)} & 201.8 & 210.5\tabularnewline
\textbf{Scaffolds:} & &\tabularnewline
\textbf{Number} & 8 188 & 2 084\tabularnewline
\textbf{Longest (Mb)} & 4.26 & 70.7\tabularnewline
\textbf{N50 (Mb)/L50} & 0.629/87 & 16.4/4\tabularnewline
\textbf{N90 (Mb)/L90} & 0.0150/1422 & 7.29/11\tabularnewline
\textbf{Average length} & 30 880 & 107 864\tabularnewline
\textbf{Number \textgreater{}10 kb} & 1 928 & 213\tabularnewline
\textbf{\% of genome \textgreater{}10 kb} & 92.5 & 97.7\tabularnewline
\textbf{Gene models} & 13 088 & 13 601\tabularnewline
\textbf{Transcripts} & 13 587 & 16 299\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2: Canadian MPB average nucleotide diversity,
Tajima's\emph{D} and \emph{F\textsubscript{ST}} estimates.}

We used Popoolation2 (Kofler et al., 2011) and the final female genome
assembly to evaluate nucleotide diversity (\selectlanguage{greek}π\selectlanguage{english}), differentiation
(\emph{F\textsubscript{ST}} ) and Tajima's~\emph{D} between northern and
southern pool-seq datasets\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
\textbf{Partition}{} & \selectlanguage{greek}\textbf{Diversity (π)North}\selectlanguage{english}{} & \selectlanguage{greek}\textbf{Diversity
(π)South}\selectlanguage{english}{} & \textbf{Tajima's~\emph{D}North}{} &
\textbf{Tajima's~\emph{D}South}{} &
\textbf{\emph{F}}\textsubscript{\textbf{\emph{ST}}}\textbf{(North vs
South)}{}\tabularnewline
Ancestral-autosomal region of neo-X{} & 1.267\%{} & 1.265\%{} & -0.965{}
& -1.050{} & 0.038{}\tabularnewline
Ancestral-X region of neo-X{} & 0.431\%{} & 0.569\%{} & -1.214{} &
-0.913{} & 0.096{}\tabularnewline
\textbf{Autosomes}{} & 1.016\%{} & 1.194\%{} & -0.790{} & -0.813{} &
0.067{}\tabularnewline
~~~~Chromosome 2{} & 0.741\%{} & 0.859\%{} & -0.880{} & -0.950{} &
0.063{}\tabularnewline
~~~~Chromosome 3{} & 0.996\%{} & 1.195\%{} & -0.793{} & -0.805{} &
0.068{}\tabularnewline
~~~~Chromosome 4{} & 1.301\%{} & 1.587\%{} & -0.599{} & -0.496{} &
0.092{}\tabularnewline
~~~~Chromosome 5{} & 0.986\%{} & 1.197\%{} & -0.824{} & -0.838{} &
0.068{}\tabularnewline
~~~~Chromosome 6{} & 0.926\%{} & 1.117\%{} & -0.791{} & -0.895{} &
0.064{}\tabularnewline
~~~~Chromosome 7{} & 0.949\%{} & 1.096\%{} & -0.738{} & -0.784{} &
0.062{}\tabularnewline
~~~~Chromosome 8{} & 0.838\%{} & 1.021\%{} & -0.804{} & -0.789{} &
0.068{}\tabularnewline
~~~~Chromosome 9{} & 1.077\%{} & 1.203\%{} & -0.846{} & -0.916{} &
0.055{}\tabularnewline
~~~~Chromosome 10{} & 1.191\%{} & 1.431\%{} & -0.867{} & -0.879{} &
0.059{}\tabularnewline
~~~~Chromosome 11{} & 1.256\%{} & 1.547\%{} & -0.840{} & -0.873{} &
0.062{}\tabularnewline
~~~~Chromosome 12{} & 1.723\%{} & 2.079\%{} & -0.749{} & -0.868{} &
0.058{}\tabularnewline
\bottomrule
\end{longtable}



\textbf{Figure Legends}

\textbf{Figure 1: Genome assembly pipeline}

Each major step en route to the final annotated assemblies is shown in
blue text. Sources of the input data (NCBI BioSample, assembly, and SRA
accessions) are indicated in parentheses. Dashed lines between the two
assembly paths indicate sharing of scaffolding information only; solid
lines indicate sharing of scaffolding and sequencing information. For
annotation of the genomes with Maker, empirical sequence sources and
the\emph{ab initio} gene predictors used are indicated.

\textbf{Figure 2: Example of manual scaffold editing based on linkage
map information}

Example of manual repositioning of regions of the scaffolds from the
Hi-C step using information from the linkage map and male versus female
dot-plots. (A) ALLMAPS graphs for a male Hi-C scaffold showing position
of markers derived from male (orange) and female (green) mapping
population samples. Sections with crossing marker positions (left
subgraph) and negative slope (right subgraph) indicate regions in the
scaffold with orientation inconsistent with the linkage map. Sections in
right subgraph that are parallel but not collinear with the trend of the
curve indicate regions in the scaffold with ordering inconsistent with
the linkage map; (B) ALLMAPS graphs for a female Hi-C scaffold showing
position of markers derived from male and female mapping population
samples; (C) Dot-plot (mummerplot with data generated by NUCmer) of the
male and female Hi-C scaffolds. Blue arrows indicate the position and
orientation of the Chicago scaffolds that were inconsistent with linkage
map information. Black arrows indicate positional reordering of these
scaffolds, and red arrows indicate a flip of the scaffold; (D) ALLMAPS
graphs for the same male Hi-C scaffold after manual repositioning
showing consistency with the linkage map information; (E) ALLMAPS graphs
for the same female Hi-C scaffold after manual repositioning, showing
consistency with the linkage map information; (F) Dot-plot of the male
and female Hi-C scaffolds after manual repositioning, showing
consistency between sexes. Blue arrows demark the position and
orientation of the Chicago HiRise scaffolds after the manual
repositioning.

\textbf{Figure 3: Linkage groups representing putative chromosomes for
the female-aligned (A) and male-aligned (B) linkage maps}

The vertical axis on the left-hand side of the figure indicates the cM
size of each linkage group, and the black horizontal lines in each
linkage group depict the positions of SNPs used to construct the maps.
The background colour of each linkage group shows recombination density;
regions of higher density (exhibiting greater linkage) are indicated by
warm colours, and regions with lower density are represented by cooler
colours. Note that the female-aligned linkage groups 10 and 11 are
labelled in reverse order from the equivalent male-aligned linkage
groups, indicated by the arrows. Chromosomes in the genome assemblies
are numbered from largest to smallest, and the assembled size of
chromosomes 10 and 11 differ between the sexes.

\textbf{Figure 4: Patterns of genome-wide differentiation and
Tajima's\emph{D} of the female mountain pine beetle genome}

(a) Manhattan plot for \emph{F\textsubscript{ST}} was calculated in 10
kb non-overlapping sliding windows with each point representing a single
estimate of \emph{F\textsubscript{ST}} ordered at the relative genomic
position and colored according to chromosome number. (b) Line graph of
Tajima's \emph{D} calculated in 10 kb non-overlapping sliding windows
for northern (red) and southern (blue) pool-seq populations of mountain
pine beetle. (c) Manhattan plot for \emph{F\textsubscript{ST}} and line
graph of Tajima's \emph{D} for neo-X. (d) Manhattan plot
for\emph{F\textsubscript{ST}} and line graph of Tajima's \emph{D} for
chromosome 4. Regions of elevated \emph{F\textsubscript{ST}} and
Tajima's \emph{D} highlighted in grey.

\textbf{Figure 5: A schematic of the male neo-X chromosome}

(A) Schematic of the neo-sex chromosomes. The neo-X chromosome
originates from the fusion of an ancestral autosome and the ancestral-X
chromosome (shown in orange). The neo-Y chromosome (shown in cyan)
originates from the other copy of the ancestral autosome found in the
neo-X. (B) Schematic of the male neo-X chromosome. Regions of cyan
indicate neo-Y scaffolds identified from the draft male assembly by
Dowle et al. (2017) that mapped on to the male neo-X chromosome in the
new genome assembly. The inset image at the top of the figure
illustrates that these scaffolds are short and not contiguous along the
chromosome, but rather interspersed among neo-X regions. The grey region
on the far-right portion of the chromosome indicates a large region
where putative neo-Y scaffolds failed to map, and likely represents the
ancestral-X portion of the chromosome.

\textbf{Supplemental Information:}

Supplemental Tables 1-10: Supplemental\_Tables\_1-10.xlsx

Supplemental Figures 1-6: Supplemental\_Figures\_1-6.pdf

Supplemental File 1: Supp\_File\_1.xlsx\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-1/Fig-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-2/Fig-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-3/Fig-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-4/Fig-4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-5/Fig-5}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

