Deforestation patterns shape population structure of the Magellanic Woodpecker (Campephilus magellanicus) in southern Chile

One important landscape-scale consequence of deforestation is reduced connectivity, which has the potential to isolate populations in ways that affect genetic diversity and population structure. Among the many regions of the world where this scenario has played out is the South American temperate forest (SATF) in southern Chile, and there is now strong concern about the population viability of forest taxa. We studied one such species, the Magellanic Woodpecker (Campephilus magellanicus), a forest specialist that is now listed as vulnerable in parts of its range in Chile. We characterized genetic variability and population structure from 33 samples of Magellanic Woodpeckers from two large but isolated populations in Nahuelbuta National Park in the Coastal mountain range and the Conguillío National Park in the Andes using ddRAD-seq method. We found lower genetic variability in Nahuelbuta than Conguillío, but inbreeding values (FIS) did not show evidence for inbreeding depression. Results suggest the presence of two genetic clusters, with an average FST value of 0.04. Phylogenetic analysis showed that the Nahuelbuta population forms a clade that is nested within the individuals from Conguillío, suggesting limited gene flow between these populations. Our results support the idea that extensive deforestation has played a role in shaping the genetic patterns that we have identified. Because of this, we emphasize the need for regional planning to increase the structural connectivity between fragments of mature native forests, to provide an opportunity for the persistence of Magellanic Woodpeckers in this region. Les modèles de déforestation façonnent la structure de la population du Pic de Magellan (Campephilus magellanicus) dans le sud du Chili RÉSUMÉ. Une conséquence importante de la déforestation à l'échelle du paysage est la diminution de la connectivité, qui peut à son tour entraîner l'isolement des populations jusqu'à affecter leur diversité génétique et leur structure. La Forêt tempérée d'Amérique du Sud dans le sud du Chili figure parmi les nombreuses régions du monde où ce scénario s'est produit, et on s'inquiète désormais fortement de la viabilité des populations des taxons forestiers qui s'y trouvent. Nous avons étudié une de ces espèces, le Pic de Magellan (Campephilus magellanicus), spécialiste des forêts maintenant classé comme vulnérable dans certaines parties de son aire de répartition au Chili. Nous avons caractérisé la variabilité génétique et la structure de la population à partir de 33 échantillons de Pics de Magellan provenant de deux populations importantes mais isolées, dans le parc national de Nahuelbuta dans la chaîne de montagnes côtières et dans le parc national de Conguillío dans les Andes au moyen de la méthode ddRAD-seq. Nous avons trouvé une variabilité génétique plus faible à Nahuelbuta qu'à Conguillío, mais les valeurs de consanguinité (FIS) n'ont pas montré de signes de dépression de la consanguinité. Nos résultats indiquent la présence de deux groupes génétiques et une valeur FIS moyenne de 0,04. L'analyse phylogénétique a révélé que la population de Nahuelbuta forme un clade imbriqué dans les individus de Conguillío, laissant entendre qu'il existe un flux génétique limité entre ces populations. Nos résultats corroborent l'idée voulant que la déforestation extensive a joué un rôle dans la formation des modèles génétiques que nous avons identifiés. Nous croyons donc qu'il est nécessaire de planifier régionalement pour accroître la connectivité structurelle entre les fragments de forêts indigènes matures, afin d'assurer la pérennité du Pic de Magellan dans cette région.


INTRODUCTION
Landscape composition and configuration can profoundly affect genetic variability (DiLeo and Wagner 2016), with serious implications for persistence, ability to adapt to changing environments, and ultimately evolutionary potential of wild populations (Fisher 1958). When populations are sufficiently large and connected, high levels of gene flow promote high genetic variability and homogenization (England et al. 2003). However, in fragmented landscapes, populations can become isolated and genetic variability may decline to the point of increasing risk of local extinction (Lande 1994, Frankham 1995, Frankham et al. 2017a) and/or promoting differentiation and the emergence of new lineages (Wright 1978). The negative consequences of isolation seem most acute for forest ecosystems, possibly because most of the species that inhabit forests are not be able to cross deforested areas (e.g., Turner and Corlett 1996, Dale 2001, Foley et al. 2005, Bruggeman et al. 2010, Frankham et al. 2017b, Botero-Delgadillo et al. 2020. Therefore, information about genetic variability of forest-associated species may be crucial to identify meaningful biological units (e.g., subspecies), estimate population viability, and ultimately, to develop conservation strategies (Walters 1991, Allendorf and Luikart 2007, Palsbøll et al. 2007. Within forest ecosystems, woodpeckers are known to be especially vulnerable to deforestation and forest degradation because of their high degree of specialization on old growth trees, which they use both to search for food and to build their nests (Lammertink 2014). Research shows that the global species richness of woodpeckers is strongly linked with tree cover, precipitation (Ilsøe et al. 2017), and presence of deadwood at the early decay stage (Tremblay et al. 2009); therefore, events that modify these conditions, such as climate change and human activities, have a strong impact on the woodpeckers' persistence (Billerman et al. 2016). Important examples of this are the cases of the two largest Campephilus species, the Ivory-billed Woodpecker (C. principalis Linnaeus 1758) and the Imperial Woodpecker (C. imperialis Gould 1832), both old-forest specialists driven to extinction as a result of habitat loss (Lammertink andEstrada 1995, Fitzpatrick et al. 2006). Likewise, research on other woodpeckers shows that habitat loss and fragmentation reduced genetic diversity, increased inbreeding, and lowered effective population sizes, e.g., Red-cockaded Woodpecker (Picoides borealis; Reed et al. 1988, Stangel et al. 1992, Haig et al. 1993, Blackwell et al. 1995, Schiegg et al. 2006, Bruggeman et al. 2010, White-backed Woodpecker (Dendrocopos leucotos; Ellegren et al. 1999), and Wryneck (Jynx torquilla; Mermod et al. 2009;Assandri et al. 2018).
The Magellanic Woodpecker (Campephilus magellanicus), the largest extant member of Campephilus (40 cm, 276-363 g; Short 1982) is an endemic and resident species of the South American temperate forests (SATF; Vuilleumier 1985). As a forest specialist, the Magellanic Woodpecker excavates mature and decayed trees both to forage for saproxylic insect larvae  and to build cavities for nesting and refuge (Ojeda 2004). Although this species is the largest of the four primary excavators of the SATF (Colaptes pitius, Veniliornis lignarius, and Pygarrhychas albogularis) that help support a diverse community of secondary cavity nesters (Altamirano et al. 2017), it might also serve as a potential seed disperser for non-Nothofagaceae trees in the Valdivian ecoregion (Soto et al. 2018). For these reasons, the Magellanic Woodpecker has been proposed as a keystone and indicator species for the SATF (Ojeda and Chazarreta 2014) and also a charismatic species because of its great social-cultural value (Arango et al. 2007).
Unfortunately, a century of landscape change and deforestation have threatened the viability of Magellanic Woodpecker populations, thus prompting interest in its ecology (Ojeda 2004, Ojeda and Chazarreta 2006, 2017, 2019. In addition to shrinking distribution due to deforestation (Ojeda andChazarreta 2014, Vergara et al. 2017), populations suffer low reproductive output (Chazarreta et al. 2012), intraspecific competition , and harm from invasive species (Jiménez et al. 2014). Presently, the species is listed as endangered and vulnerable on the northern and southern extents of its distribution in Chile, respectively (Servicio Agrícola y Ganadero 2015). In this context, a better understanding of the genetic diversity and population structure of the Magellanic Woodpecker is crucial to develop conservation plans aimed to maintain genetic variation, particularly in small and/or isolated populations, and also to prevent the fixation of deleterious alleles with potentially negative fitness consequences (Van Dyke 2008).
The most isolated population of Magellanic Woodpeckers (> 100 km apart from other local populations) is located in the northern extent of their distribution, within the boundaries of the Nahuelbuta National Park (hereafter Nahuelbuta; 6832 ha; 37°4 9′ 03″ S, 73° 02′ 03″ W; Vergara et al. 2017). This population is located within the northern extent of the Valdivian rainforest ecoregion, one of the world's 32 biodiversity hotspots (Myers et al. 2007). This hotspot has been highly threatened by historic deforestation in the last centuries on the lowlands and low-slope hills, leading to a discontinuity of the forest between the two mountain ranges dominating the south-central part of Chile ( Fig.1; Olson et al. 2001, Otero 2006. Consequently, Nahuelbuta is now an island of relatively continuous old-growth forest within a heavily fragmented region in the Nahuelbuta mountain range; extensive stretches of unsuitable habitat separate Nahuelbuta from mature forests in the Andes ( Fig.1; Dinerstein et al. 1995, Cox andKesley 2012).
The closest known potential nesting habitat for Magellanic Woodpecker to Nahuelbuta is located in high slope forests in the proximity of the Malleco National Reserve (38° 08′ 03″ S, 71° 47′ 05″ W), about 100 km to the east on the Andes mountain range. However, the closest relatively large local population inhabiting low-slope mature forests can be found ~50 km south from Malleco National Reserve, within the boundaries of the Conguillío National Park (hereafter Conguillío; 4500 ha; 38° 38′ 39″ S, 71°3 8′ 42″ W; Fig. 1). Conguillío features similar tree species composition when compared to Nahuelbuta, with mixed and pure Nothofagaceae and Araucaria araucana, but is relatively well connected with other similar forests along the Andes (Aguilera-Betti et al. 2017, Soto 2019).
To help inform the conservation status of this species, we compared samples of Magellanic Woodpeckers from Nahuelbuta and Conguillío national parks, using measures of genetic variability, inbreeding, and population structure, to understand the possible genetic consequences of historic deforestation in this region. We predicted that (1) isolated populations in Nahuelbuta distribution. The left map shows the two sampled national parks and the remnant native forests, illustrating the extensive deforestation in the lowlands between the two main mountain ranges. Other land-cover types such as agriculture, urban areas, and exotic plantations are labeled as 'Other'. The upper-right map shows the Magellanic Woodpecker's distribution in southern South America (grey area) and the northern extent of the range where both national parks are located. The lower-right map shows the Chilean remnants of native forests in the northern extent of the South American temperate forests and the geographical extent of the left panel.
would have lower genetic diversity and greater levels of inbreeding than Conguillío, and (2) our two populations would show signs of genetic divergence.

Study sites and sampling
A total of 33 Magellanic Woodpecker individuals were sampled during austral spring (September to December) of 2017 in Nahuelbuta (N = 14) and Conguillío (N = 19), central-southern Chile (Fig. 1). Nahuelbuta is the only protected area within the coastal mountain range that includes A. araucana, covering 61.91 km² (98.2%) of native forest and includes shrub vegetation in degraded forest. Conguillío (545 km²), in the Andes Mountains, presents a native forest cover of 230.72 km² (42.3%) and other species such as cushion plants, mosses, perennial grasses, scrub, and lichens growing on volcanic rocks (CONAF 2019).
Birds were captured using mist nets (18 m long × 3 m high) and marked with uniquely numbered metal bands (National Band and Tag Co., Newport, Kentucky, USA; Porzana Ltd, UK, provided by Servicio Agrícola y Ganadero de Chile). Before release, a small blood sample was obtained by brachial venipuncture (ca. 75 µL) and stored on FTA cards (Whatman®) for subsequent molecular analysis.

DNA sequencing
We generated single nucleotide polymorphisms (SNPs) through a double-digest restriction site-associated DNA sequencing (ddRAD-seq) method, following the protocol described in Thrasher et al. (2018). We isolated genomic DNA with the DNeasy Blood and Tissue kit (Qiagen, Valencia, CA, USA) following the manufacturer's protocols. Briefly, we digested ~200 ng genomic DNA from each individual with the enzymes SbfI and MspI, and ligated adapters on both the 5′ and 3′ ends. After digestion and ligation, these samples were pooled in groups of 19-20, based on unique barcodes present in the 5′ adapters. We subsequently selected DNA fragments in the 400-700 bp range for each pool separately using BluePippin (Sage Science). We added Illumina index groups and enriched our libraries through PCR. Finally, we combined each pool in equimolar ratios to create a single library for sequencing (together with samples from other projects) on one lane of an Illumina NextSeq 500 in mid output mode (single-end 151 bp). Sequencing was conducted at the Cornell Institute of Biotechnology.

Bioinformatics and assembly of RAD loci
We obtained an average of 0.57 ± 0.2 million reads per individual.
We assessed read quality using FastQC version 0.11.6 (https:// www.bioinformatics.babraham.ac.uk/projects/fastqc/) and decided to exclude the last 4 bp on the 3′ end, which tended to show comparatively lower average quality calls by trimming sequences to 147 bp using fastX trimmer (Gordon and Hannon 2017). We subsequently discarded reads with bases below a Phred quality score of 10 (90% call accuracy) or with more than 5% between 10 and 20 (99% call accuracy) using fastq_quality_filter (fastx-Toolkit). We demultiplexed the reads using the process_radtags module from the stacks pipeline version 1.48 (Catchen et al. 2013) obtaining files with sequences exclusive to each individual. All sequences were trimmed to 140 bp, the length of the shortest sequences after the 5-7 bp inline barcodes used for multiplexing was removed. We assembled sequences into RAD loci using the de novo pipeline from stacks (ustacks/cstacks/sstacks), requiring a minimum depth of coverage of 5 reads (m parameter) and allowing up to 5 mismatches between different alleles within individuals (M parameter) or among individuals (n parameter; Thrasher et al. 2018). We ran the rxstacks error correction module and a final iteration of cstacks/sstacks. This bioinformatics pipeline generated a catalog of 21,716 RAD loci from which we exported SNPs in different formats for downstream analyses using the population module from stacks, allowing up to 20% missing data per locus and requiring a minimum depth of coverage of 10 ×. We exported both all the SNPs from the RAD loci that passed our filtering criteria (2669 SNPs) and only the first SNP from each rad locus to avoid including tightly linked markers (1674 SNPs).

Population genomic analyses
We calculated summary statistics for each population and FST between populations using the population module in Stacks. We used four methods to analyze our SNP data: a principal component analysis (PCA) with the SNPrelate R package version 3.3 (Zheng et al. 2012), the clustering software STRUCTURE version 2.3.4 (Pritchard et al. 2000), the haplotype-based clustering analysis implemented in fineRADstructure v 0.2 (Malinsky et al. 2018), and we built a Maximum Likelihood phylogenetic tree in RaxML version 8.2.4 (Stamatakis 2014). For the principal component analysis and the phylogenetic tree, we used the dataset of 2669 SNPs, and for the STRUCTURE analysis, we used a single SNP per RAD locus (total of 1674 SNPs) to avoid using highly linked SNPs. The fineRADstructure analysis used haplotype information from 7207 variant and invariant RAD loci, which contained the 2669 SNPs used in other analyses. For the STRUCTURE analysis, we implemented the admixture ancestry model and used correlated allele frequencies, running the program for 750,000 generations, discarding the first 250,000 as burn-in. We conducted 10 runs exploring values of K from 1 to 5, with a different random seed per iteration. We determined the most likely value of K in STRUCTURE HARVESTER version 0.6.94 (Earl and vonHoldt 2012), following the methods described by Evanno et al. (2005). Different iterations from the same K value were combined in CLUMPP version 1.1.2 (Jakobsson and Rosenberg 2007) and displayed using DISTRUCT version 1.1 (Rosenberg 2004). We ran the fineRADstructure analysis under default parameters, excluding one individual from Conguillío with higher than average missing data. We built the Maximum Likelihood phylogenetic tree from our SNP data by implementing the ASC_GTRGAMMA model and the Lewis correction for ascertainment bias. We included 3 individuals from a distant population from Navarino Island (sequenced for a different project) as the outgroup to root the phylogeny and conducted 500 bootstrap replicates to assess node support. This analysis uses loci in which both alleles are found in homozygosity in at least 1 individual, allowing us to retain 1078 SNPs in our final alignment. Finally, we attempted to reconstruct the demographic history of the two populations using (G-PhoCS) version 1.2.2 (Gronau et al. 2011), yet despite using various strategies to conduct the analysis (e.g., subsampling individuals or excluding migration estimates), we did not observe convergence in the different parameters when inspecting the traces in Tracer 1.7 (Rambaut et al. 2018). It is possible that the divergence among these populations is too shallow for the model to accurately estimate the different demographic parameters and we, therefore, do not discuss these results.

RESULTS
We recovered 2669 SNPs in 1674 RAD loci across all individuals. Private alleles (PA), percent polymorphic loci (% PL), nucleotide diversity (Pi), and both observed (Ho) and expected (He) heterozygosity levels were lower in Nahuelbuta compared to Conguillío (Table 1). This pattern remained after randomly subsampling individuals from the Conguillío population to the same number of individuals from the Nahuelbuta population (Appendix 1). For each population, the expected and observed levels of heterozygosity were similar. Consequently, FIS values were close to zero and similar for both populations (Table 1) suggesting that these populations are not inbred. Our samples formed two clusters, coinciding with the sampling sites, in a PCA derived from 2669 SNPs. The first two principal components (PC1 and PC2) explained 7.03 % and 5.10 % of the total variance, respectively. The PC1 shows a significant separation of the two subpopulations (i.e., Nahuelbuta and Coguillío; Fig. 2A). The STRUCTURE analysis derived from 1674 SNPs (a single SNP per RAD locus) was consistent with the PCA in showing differentiation between Nahuelbuta and Conguillío, yet some individuals showed low levels of admixture (Fig. 2B). Higher values of K were not supported in our Structure analysis (Appendix 2). Results from fineRADstructure, a haplotype-based analysis, also supported the genetic differentiation between Nahuelbuta and Coguillío, showing signals of increased coancestry within each subpopulation (Appendix 3). The average FST value between Nahuelbuta and Conguillío was 0.04 (range 0-0.66 across all loci; Fig. 2C) suggesting a relatively low differentiation between both populations. Finally, phylogenetic analysis using our SNP data, shows that the Nahuelbuta population forms a clade that is nested within the individuals from Conguillío (Fig. 3), a pattern consistent with the Nahuelbuta population being derived from the Conguillío population.

DISCUSSION
The genetic variability of populations depends upon a wide range of factors, including dispersal capacity (Rappole 2013), the degree of ecological specialization of the species (Anders et al. 1998, Crochet 2000, Harris and Reed 2002, and landscape attributes (Bélisle andSt. Clair 2001, McRae 2006). In this study, we show differences in genetic variability and evidence of population structure between Magellanic Woodpeckers occupying Nahuelbuta and Conguillío national parks in the northern extent of the Valdivian ecoregion in Chile. Specifically, we found a low number of private alleles and polymorphic loci as well as reduced levels of nucleotide diversity and heterozygosity in Nahuelbuta (Table  1), which are consistent with lower genetic variability.
Although Magellanic Woodpeckers can fly long distances and overcome potential barriers to dispersal (Winkler et al. 2016), they also are reluctant to cross open areas and venture into a deforested matrix. Research on the movement behavior of Magellanic Woodpeckers shows that dominant males determine their movement trajectory and patch residence times according to observable and memorized information of local habitat quality (Vergara et al. 2015, suggesting that individuals avoid the use of a hostile matrix to perform dispersal, similarly to other woodpecker species (e.g., Pasinelli andWalters 2002, Trainor et al. 2013). Hence, evidence suggests that lower genetic variability values observed in Nahuelbuta may be a consequence of the open matrix surrounding this area, which would limit the ability of the forest-dependent woodpeckers to disperse (Gorman 2014, Lammertink 2014). Although we did not find evidence for inbreeding (FIS > 0) in these populations (Table 1), it is important to note that inbreeding is a gradual phenomenon and likely to increase over time in small and/or closed populations in which mating individuals will have, eventually, some degree of kinship (Ralls et al. 2014, Chen et al. 2016. Across time, inbreeding depression, defined as the decline in fitness observed in inbred progeny (Keller and Waller 2002), may occur in these populations because of the overall increase in homozygosity and exposure of recessive deleterious alleles (Charlesworth and Charlesworth 1987). Evidence in birds has shown that inbreeding depression negatively affects survival and reproductive success (Bensch et al. 1994, Spottiswoode and Møller 2004, Cassinello 2005, Szulkin et al. 2007, Hemmings et al. 2012. Indeed, several studies have concluded that dispersal barriers often promote inbreeding over time (Van Tienderen and Van Noordwijk 1988, Daniels and Walters 2000, Glémin et al. 2003, Szulkin and Sheldon 2008. Consequently, if Magellanic Woodpeckers avoid using a deforested and presumably hostile matrix, inbreeding may gradually rise in the Nahuelbuta population to the point of reducing fitness. The low level of genetic structure detected between Nahuelbuta and Conguillío samples suggests that divergence between both clusters is relatively recent (Freeland et al. 2012). Unfortunately, demographic models were not able to infer the divergence time between both clusters, likely because of low differentiation between populations and the long generation time of the species (Lois et al. 2020). However, the phylogenetic tree showed that the Nahuelbuta clade is nested within the individuals from Conguillío, indicating recent isolation. Deforestation during recent decades has likely reduced connectivity among remnant forests in ways that affect population structure ).
Collectively, our results suggest that Nahuelbuta woodpeckers seem to be in an initial phase of loss of genetic variability, likely due to their reluctance to cross a deforested matrix (Lammertink 2014). If Magellanic Woodpeckers remain isolated in Nahuelbuta, we predict that genetic variability will further decline and may eventually compromise their ability to adapt to environmental change.
According Within this context, our results provide further evidence that the persistence of the Magellanic Woodpeckers in the Nahuelbuta Mountain range remains uncertain. We recommend that conservation efforts focus on engaging the public and private sectors to protect and restore native mature forests in ways that increase structural connectivity at regional or national scales.  Table A1. Population genetic parameters obtained after randomly subsampling Conguillío individuals to match the sample size for Nahuelbuta (N=14). We obtain similar results to those reported in Table 1. Note that although we did not change the number of samples from Nahuelbuta, the values of some parameters changed slightly because the missing data filter is now applied to a different total number of individuals (and we only exported loci present in at least 80% of individuals  Figure A2. Exploration of the K values in Magellanic woodpecker samples. (A) Structure plots values of K from ranging from 2 to 5. (B) Likelihood and delta K sensu Evanno et al. 2005. Values of K beyond 2 are not supported. Figure A3. fineRADstructure plot showing higher levels of co-ancestry between individuals within each population than between individuals from different populations.