Code from: Novel Arctic colonization origin of the Arctic lineage in Phyllodoce (Ericaceae)
Data files
Sep 23, 2026 version files 171.17 KB
-
1_WGS_sequene_processing.R
52.14 KB
-
2_Phylogenetic_analyses.R
6.38 KB
-
3_Estimation_of_demographic_inference.R
38.57 KB
-
4_Ecological_niche_modeling.R
25.52 KB
-
5_Validation_with_Simulation_Data.R
23.19 KB
-
6_Pairwise_Sequentially_Markovian_Coalescent.R
13.29 KB
-
README.md
12.08 KB
Abstract
The Arctic flora is widely considered evolutionarily young, reflecting the emergence of unprecedentedly cold environments during the late Pliocene–early Pleistocene transition. Although many Arctic plants are thought to be derived from temperate alpine ancestors, the historical relationship between high-latitude occupancy and lineage divergence remains unclear. Specifically, it remains unclear whether Arctic lineages originated from ancestors that already occupied high-latitude regions before the Arctic environment formed or originated through colonization after the Arctic environment was established. To address this question, we test demographic models associated with the divergence of the Arctic subgenus Eu-Phyllodoce and its sister alpine subgenus (Parabryanthus) within Phyllodoce (Ericaceae). Whole-genome data from Phyllodoce were used to infer phylogenetic relationships and demographic history. Demographic parameters were estimated in an isolation-with-migration model and evaluated using log-likelihood ratio tests. Historical demographic trajectories were additionally reconstructed using the pairwise sequentially Markovian coalescent (PSMC). Demographic inference recovered a smaller ancestral effective population size for Eu-Phyllodoce than for Parabryanthus, indicating disproportional retention of ancestral polymorphism. PSMC further showed pronounced demographic expansion in Eu-Phyllodoce following divergence of ancestral lineages, consistent with a strong founder effect in the Arctic lineage. These results are consistent with the scenario that Eu-Phyllodoce originated through novel Arctic colonization. Our findings support the view that many Arctic plant lineages are evolutionarily young and highlight the roles of climatic cooling during the late Pliocene–early Pleistocene and of novel Arctic colonization in shaping the assembly of the extant Arctic flora.
Dataset DOI: to be assigned
Overview
This repository contains R scripts used to process whole-genome sequencing (WGS) data and to perform phylogenetic, demographic, and ecological niche analyses for the study:
Novel Arctic colonization origin of the Arctic lineage in Phyllodoce (Ericaceae)
All analyses were conducted using reproducible, script-based workflows. Each script is numbered to reflect the recommended order of execution.
Script File Structure
All scripts are written in R and organized as follows:
- 1_WGS_sequene_processing.R
- 2_Phylogenetic_analyses.R
- 3_Estimation_of_demographic_inference.R
- 4_Ecological_niche_modeling.R
- 5_Validation_with_Simulation_Data.R
- 6_Pairwise_Sequentially_Markovian_Coalescent.R
All required R packages are loaded within each script.
Computational Environment
Analyses were performed under the following software and system configurations:
Shells
- zsh 5.9 (x86_64-apple-darwin23.0)
- bash 5.2.37(1)-release (aarch64-apple-darwin23.4.0)
R
- R version 4.4.2 (2024-10-31)
- RStudio version 2024.09.0+375
Conda
- conda version 25.1.1
External software dependencies (e.g., RAxML, ASTRAL, MACSE, IMa3) are called from within the scripts and must be installed and accessible via the system PATH.
(Note: the RAxML-NG graphical user interface is not required.)
Script Descriptions
1. 1_WGS_sequence_processing.R
Description:
This script processes raw whole-genome sequencing (WGS) data into locus-specific FASTA files for downstream phylogenetic and demographic analyses.
Main steps:
- R environment setup
- Terminal and PATH configuration
- Quality control of raw FASTQ files
- Conversion of FASTQ files to SAM, BAM/BAI, and haplotype-resolved FASTA formats
- FASTA file processing and management
- Conversion from per-sample FASTA files to per-locus FASTA files
- Sequence alignment using MACSE (configured for Apple Silicon / arm64 systems)
- Alignment quality assessment (inspection of frameshifts, stop codons, and missing data)
- Filtering of single-copy nuclear genes
2. 2_Phylogenetic_analyses.R
Description:
This script performs phylogenetic inference using both summary coalescent and maximum-likelihood approaches.
Main steps:
- Maximum-likelihood tree inference using RAxML-HPC
- Preparation of gene trees and inference of species trees using weighted ASTRAL (wASTRAL)
- Construction of concatenated alignments for maximum-likelihood analyses
3. 3_Estimation_of_demographic_inference.R
Description:
This script estimates divergence times and historical demographic parameters using the Isolation-with-Migration framework implemented in IMa3.
Main steps:
- Preparation of input files for IMgc
- Execution of IMgc and formatting of input for IMa3
- IMa3 analyses using M-mode
- IMa3 analyses using L-mode
- Extraction and summarization of demographic estimates from M-mode
- Extraction and summarization of demographic estimates from L-mode
4. 4_Ecological_niche_modeling.R
Description:
This script performs ecological niche comparisons between Eu-Phyllodoce and Parabryanthus using the ecospat framework.
Main steps:
- Preparation of occurrence data
- Preparation of bioclimatic variables
- Extraction and screening of environmental variables
- Preparation of background environmental data
- Ecological niche visualization using PCA-based plots
- Statistical niche equivalency and similarity tests
- Visualization of niche overlap, PCA contributions, and density plots
5. 5_Validation_with_Simulation_Data.R
Description:
This script performs validation of demographic parameters estimated in empirical data using simulated sequences from msprime.
Main steps:
- Simulation of sequence files
- Check compatibility with infinite site model
- Prepare files for IMa3
- Run IMa3
- Plot results
6. 6_Pairwise_Sequentially_Markovian_Coalescent.R
Description:
This script performs Pairwise Sequentially Markovian Coalescent (PSMC) analysis of each WGS sample to estimate the historical demographic trajectory per sample using psmcr framework.
Main steps:
- Conversion of FASTQ files to SAM, BAM/BAI, vcf files
- Run PSMC
- Prepare dataframe for plots
- Plot PSMC result
- Save summary of demographic curves
Notes
- Scripts are intended to be executed sequentially (1–6).
- File paths and system-specific settings may require modification depending on the user’s local environment.
- Parameter values and analytical settings correspond to those reported in the associated manuscript.
Software Citations
The following software packages and programs were used in data processing, phylogenetic inference, demographic analyses, and ecological niche modeling.
-Aiello-Lammens M.E., Boria R.A., Radosavljevic A., Vilela B., Anderson R.P. 2015. spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography. 38:541–545.
- Andrews S. 2010. FastQC: A quality control tool for high throughput sequence data. Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
- Barrett T., Dowle M., Srinivasan A., Gorecki J., Chirico M., Hocking T., Schwendinger B., Krylov I. 2025. data.table: Extension of
data.frame[R package version 1.17.0]. Available from: https://CRAN.R-project.org/package=data.table. - Baumdicker F., Bisschop G., Goldstein D., Gower G., Ragsdale A.P., Tsambos G., Zhu S., Eldon B., Ellerman E.C., Galloway J.G., Gladstein A.L., Gorjanc G., Guo B., Jeffery B., Kretzschumar W.W., Lohse K., Matschiner M., Nelson D., Pope N.S., Quinto-Cortés C.D., Rodrigues M.F., Saunack K., Sellinger T., Thornton K., van Kemenade H., Wohns A.W., Wong Y., Gravel S., Kern A.D., Koskela J., Ralph P.L., Kelleher J. 2022. Efficient ancestry and mutation simulation with msprime 1.0. Genetics. 220:iyab229.
- Bougeard S., Dray S. 2018. Supervised Multiblock Analysis in R with the ade4 Package. J. Stat. Softw. 86:1–17.
- Chamberlain S., Barve V., Mcglinn D., Oldoni D., Desmet P., Geffert L., Ram K. 2025. rgbif: Interface to the Global Biodiversity Information Facility API [R package version 3.8.1]. Available from: https://CRAN.R-project.org/package=rgbif.
- Chamberlain S.A., Boettiger C. 2017. R Python, and Ruby clients for GBIF species occurrence data. PeerJ PrePr. 5:e3304v1.
- Danecek P., Bonfield J.K., Liddle J., Marshall J., Ohan V., Pollard M.O., Whitwham A., Keane T., McCarthy S.A., Davies R.M., Li H. 2021. Twelve years of SAMtools and BCFtools. GigaScience. 10:giab008.
- Di Cola V., Broennimann O., Petitpierre B., Breiner F.T., D’Amen M., Randin C., Engler R., Pottier J., Pio D., Dubuis A., Pellissier L., Mateo R.G., Hordijk W., Salamin N., Guisan A. 2017. ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography. 40:774–787.
- Dray S., Dufour A.-B. 2007. The ade4 Package: Implementing the Duality Diagram for Ecologists. J. Stat. Softw. 22:1–20.
- Fick S.E., Hijmans R.J. 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37:4302–4315.
- Hester J., Wickham H., Csárdi G. 2025. fs: Cross-Platform File System Operations Based on “libuv” [R package version 1.6.6]. Available from: https://CRAN.R-project.org/package=fs.
- Hey J., Wang K. 2019. The effect of undetected recombination on genealogy sampling and inference under an isolation-with-migration model. Mol. Ecol. Resour. 19:1593–1609.
- Hijmans R.J. 2025a. terra: Spatial Data Analysis [R package version 1.8-42]. Available from: https://CRAN.R-project.org/package=terra.
- Hijmans R.J. 2025b. raster: Geographic Data Analysis and Modeling [R package version 3.6-32]. Available from: https://CRAN.R-project.org/package=raster.
- Hijmans R.J., Phillips S., Leathwick J., Elith J. 2024. dismo: Species Distribution Modeling [R package version 1.3-16]. Available from: https://CRAN.R-project.org/package=dismo.
- Langmead B., Salzberg S.L. 2012. Fast gapped-read alignment with Bowtie 2. Nat. Methods. 9:357–359.
- Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R., 1000 Genome Project Data Processing Subgroup. 2009. The Sequence Alignment/Map format and SAMtools. Bioinform. 25:2078–2079.
- Massicotte P., South A. 2023. rnaturalearth: World Map Data from Natural Earth [R package version 1.0.1]. Available from: https://CRAN.R-project.org/package=rnaturalearth.
- Morgan M., Ramos M. 2024. BiocManager: Access the Bioconductor Project Package Repository [R package version 1.30.25]. Available from: https://CRAN.R-project.org/package=BiocManager.
- Müller K., Wickham H. 2023. tibble: Simple Data Frames [R package version 3.2.1]. Available from: https://CRAN.R-project.org/package=tibble.
- Pagès H., Aboyoun P., Gentleman R., DebRoy S. 2024. Biostrings: Efficient manipulation of biological strings [R package version 2.74.1]. Available from: https://bioconductor.org/packages/Biostrings.
- Paradis E. 2010. pegas: an R package for population genetics with an integrated–modular approach. Bioinform. 26:419–420.
- Paradis E, Li H. 2023. psmcr: Pairwise sequentially markovian coalescent [R package version 0.1-4]. Available from: https://github.com/emmanuelparadis/psmcr.
- Pfeifer B., Wittelsbürger U., Ramos-Onsins S.E., Lercher M.J. 2014. PopGenome: an efficient Swiss army knife for population genomic analyses in R. Mol. Biol. Evol. 31:1929–1936.
- Posit team. 2024. RStudio: Integrated Development Environment for R. Available from: http://www.posit.co/.
- R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.
- Ranwez V., Douzery E.J.P., Cambon C., Chantret N., Delsuc F. 2018. MACSE v2: Toolkit for the Alignment of Coding Sequences Accounting for Frameshifts and Stop Codons. Mol. Biol. Evol. 35:2582–2584.
- SAMtools. 2020. BCFtools. Available from: https://github.com/samtools/bcftools/ 1.11.
- South A., Michael S., Massicotte P. 2024. rnaturalearthdata: World Vector Map Data from Natural Earth Used in “rnaturalearth” [R package version 1.0.0]. Available from: https://CRAN.R-project.org/package=rnaturalearthdata.
- Stamatakis A. 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinform. 30:1312–1313.
- Wickham H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York.
- Wickham H. 2023. stringr: Simple, Consistent Wrappers for Common String Operations [R package version 1.5.1]. Available from: https://CRAN.R-project.org/package=stringr.
- Wickham H., François R., Henry L., Müller K., Vaughan D. 2023. dplyr: A Grammar of Data Manipulation [R package version 1.1.4]. Available from: https://CRAN.R-project.org/package=dplyr.
- Woerner A.E., Cox M.P., Hammer M.F. 2007. Recombination-filtered genomic datasets by information maximization. Bioinform. 23:1851–1853.
- Zhang C., Mirarab S. 2022. Weighting by Gene Tree Uncertainty Improves Accuracy of Quartet-based Species Trees. Mol. Biol. Evol. 39:msac215.
