应用环境关联分析检测本地适应位点
Detecting Loci Underlying Local Adaptation Using Environmental Association Analysis   

研究背景

随着全球气候和环境的剧烈变化，理解生物适应环境变化的遗传响应，已成为进化生物学研究中的一个重要问题。利用RAD-seq、全基因组重测序等技术获得SNP (single nucleotide polymorphism, 单核苷酸多态性)，以SNP作为遗传标记进行选择信号分析，无论对于模式物种还是非模式物种，均已成为目前适应性进化遗传机制研究的主流方法(Pardo-Diaz et al., 2015)。基于全基因组范围内SNP数据，在野生群体的适应机制研究中检测受选择信号主要有两类方法：F_ST离群值检测 (outlier analysis, OA) 和环境关联分析 (environmental association analysis, EAA) (Liberles et al., 2020)，其中F_ST离群值检测的具体流程已在本手册另一篇文章 (金铃和郭宝成，2021) 中详细介绍。相较于F_ST离群值检测，环境关联分析除需基因型数据外，还需环境变量数据，通过特定的方法分析某个环境变量与SNP等位基因频率变异间的关联，从而 推断选择压力作用于哪些位点或基因组区域 (Rellstab et al., 2015)。
        环境关联分析常用的程序或软件包括BAYENV、BAYENV2、LFMM、BAYESCENV、VEGAN等，常见方法及软件的基本原理、适用的数据类型等参见附表1 (译自Rellstab et al., 2015 Table 1)。其中BAYENV及BAYENV2使用频率较高 (Ahrens et al., 2018)。BAYENV是一种贝叶斯方法，它通过指定一组中性位点估算群体间等位基因频率变异的零模型，而后检验加入环境变量的新模型相对于原模型是否可以更好地拟合数据，从而推断与环境变量相关联的遗传变异 (Coop et al., 2010; Guenther and Coop, 2013)。相较于其他可进行同类分析的方法，BAYENV具有可对群体间样本大小的差异和群体间共享的群体历史进行校正的优势 (Guenther and Coop, 2013)。本文将以BAYENV2软件为例，向读者介绍基于重测序数据利用BAYENV2进行环境关联分析的详细分析方案，该分析方案主要参考Guo et al. (2015, 2016a) 两篇已发表的工作。该分析方案已在鱼类野生群体的适应性机制研究中广泛应用 (如Magalhaes et al., 2021)，也可应用于鱼类以外的其他生物类群 (如Guo et al., 2016b; Fraik et al., 2020)。

仪器设备

1. 服务器 (型号：Inspur NF5280M5 操作系统：CentOS Linux release 7.5.1804 (Core)；CPU：Intel(R) Xeon(R) Gold 6230 CPU @ 2.10 GHz；80核，512 G内存)
2. 普通个人电脑 (需要已安装好R)
   

软件版本信息及下载地址

1. BWA-MEM v0.7.17-r1188: https://github.com/lh3/bwa
2. SAMtools v1.8: https://github.com/samtools/samtools/releases/tag/1.8
3. BCFtools v1.8: https://github.com/samtools/bcftools/releases/tag/1.8
4. VCFtools v0.1.13: https://github.com/vcftools/vcftools/releases/tag/v0.1.13
5. Java v1.8.0_151: https://java.com/zh-CN/
6. PGDSpider 2.1.1.5: http://cmpg.unibe.ch/software/PGDSpider/
7. BAYENV2: https://bitbucket.org/tguenther/bayenv2_public/src/master/
8. R v3.6.3: https://www.r-project.org/
9. ggplot2 (R package, v3.3.3): https://cran.r-project.org/web/packages/ggplot2/index.html
10. scales (R package, v1.1.1): https://CRAN.R-project.org/package=scales
    

实验步骤

1. 产生一个高质量的SNP数据集
   利用BWA、SAMtools、BCFtools、VCFtools等软件，将二代测序平台产生的对研究群体的重测序数据比对至参考基因组上，并进行SNP检测、SNP过滤等，获得一个高质量的SNP数据集 (一般为VCF格式)。详细操作流程参见本手册另一篇文章 (金铃和郭宝成，2021)。
   
2. 格式转换
   使用PGDSpider软件将SNP数据集由VCF格式转换为BAYENV格式。
   $ java -jar PGDSpider2-cli.jar -inputfile test.vcf -inputformat VCF -outputfile test.bayenv -outputformat BAYENV -spid VCF_BAYENV.spid
   注：上述命令不加-spid参数运行，可自动生成配置文件模板template_ VCF_BAYENV.spid，需依据实际数据情况修改该模板中的内容，而后加-spid参数再次运行上述命令。
   
3. 准备环境变量文件
   在环境变量文件中，依据SNP数据集中的群体顺序，列出参与环境关联分析的环境变量的值。每个环境变量必须进行标准化，即减去平均值后除以参与分析群体的标准差。每个环境变量占一行，数值之间以Tab键分隔。具体文件格式参考BAYENV2软件提供的示例文件PCs.env (https://bitbucket.org/tguenther/bayenv2_public/src/master/example_input.zip)。
   注：常见环境变量一般可从公共数据库和网站等获取，例如可获取气候变量的常用公共资源包括IRI/LDEO Climate Data Library (http://iridl.ldeo.columbia.edu/)、WorldClim Database (https://www.worldclim.org/)等。
   
4. 估算协方差矩阵
   输入一组中性位点，估算群体间的协方差矩阵，即估算群体间等位基因频率变化的零模型，用于下一步检测环境相关性中排除掉等位基因频率变异与群体结构、环境或表型变异间假的关联。一般可将BayeScan或其他软件和方法鉴定到的受选择位点 (如F_ST离群SNPs) 从SNP数据集中剔除，再进行SNP数据集缺失数据、连锁不平衡的过滤。过滤后最终保留几千个SNPs，协方差矩阵的估算即可收敛。
          协方差矩阵的估算通常需要用不同的随机种子 (random seed) 独立运行三次以上，以确保矩阵估算已经收敛。在下一步环境关联分析中，可使用多次协方差矩阵估算结果的平均值。
   $ java -jar PGDSpider2-cli.jar -inputfile neutral_loci.vcf -inputformat VCF -outputfile neutral_loci.bayenv -outputformat BAYENV -spid VCF_BAYENV.spid
   $ bayenv2 -i neutral_loci.bayenv -p 17 -k 50000 -r 3240 > neutral_loci. matrix
   
5. 检测环境相关性
   使用BAYENV2软件自带的脚本calc_bf.sh (https://bitbucket.org/tguenther/bayenv2_public/src/master/calc_bf.sh)，估算每个SNP对每个环境变量的贝叶斯因子 (Bayes factor, BF)。输出文件中贝叶斯因子的大小，可代表环境变量与等位基因频率变异之间相关的支持度；贝叶斯因子的值越大，说明该SNP与特定环境变量之间的关联越强。
   检测环境相关性同样需要用不同的随机种子独立运行三次及以上，取多次运行结果的平均值，或多次分析均鉴定到的环境关联SNPs。
   $ calc_bf.sh test.bayenv ENVIRONFILE neutral_loci.matrix.mean 17 50000 2
   注：calc_bf.sh后的参数分别为BAYENV格式输入文件、环境变量文件、协方差矩阵文件、群体数目、MCMC链迭代数目、环境变量数目。为保证脚本正常运行，上述参数应严格按照以上顺序添加至脚本后。
   
6. 结果分析及绘图
   生成文件bf_environ.ENVIRONFILE第一列为SNP ID，后续几列为该SNP与各个环境变量间的贝叶斯因子。该文件可用于绘制SNP等位基因频率变异与环境变量间相关联的曼哈顿图或频率分布直方图。下文中图1和图2是用Guo et al. (2016a) 中17个波罗的海大西洋鲱鱼群体对盐度的环境关联分析的结果生成的 (https://datadryad.org/stash/dataset/doi:10.5061/dryad.d8b93)，供读者参考。
   1) 绘制曼哈顿图，以各个SNP ID为横坐标，以贝叶斯因子以10为底的对数为纵坐标，以log₁₀[Bayes factor (BF)] > 1.5作为判定SNP是否与该环境变量强相关的阈值 (Jeffreys, 1961) (图1)。
   加载R包ggplot2：
   > library(ggplot2)
   读取BAYENV2输出文件：
   > bf <- read.table("bf_environ.ENVIRONFILE")
   > colnames(bf) <- c("SNP_ID","Salinity")
   生成一列SNP序号：
   > bf$NUM <- c(1:nrow(bf))
   绘图：
   > man_plot <- ggplot(bf,aes(x=NUM,y=log10(Salinity)))+
     geom_point(size=2)+
     geom_hline(aes(yintercept=1.5),size=1.5,linetype="dashed")+
     xlab("SNPs")+
     ylab(expression(log[10](BF)))+
     ggtitle("Salinity")+
     theme_bw()+
     theme(
       axis.text=element_text(size=30),
       axis.title=element_text(size=35),
       plot.title=element_text(size=35,hjust=0.5),
       panel.border=element_rect(size=1.5),
       panel.grid=element_blank()
   )
   保存图像：
   > ggsave("man_plot.pdf",plot=man_plot,width=28,height=6)
   
   
   图1. SNP等位基因频率变异与环境变量相关联的曼哈顿图 虚线表示log₁₀[Bayes factor (BF)] = 1.5的阈值。
   
   2) 绘制SNP等位基因频率变异与环境变量相关联的频率分布直方图，从该图可以看出贝叶斯因子值的分布 (图2)。
   加载R包scales：
   > library(scales)
   绘图：
   > freq_plot <- ggplot(bf,aes(x=log10(Salinity)))+
     geom_histogram(bins=100,aes(y=after_stat(count/sum(count)*100)),
                   color="#BABABA",fill="#BABABA")+
     xlab(expression(log[10](BF)))+
     ylab("Frequency (percentage)")+
     ggtitle("Salinity")+
     theme_classic()+
     scale_x_continuous(limits=c(-1,1.5),oob=squish)+
     theme(axis.text=element_text(size=20),
           axis.title=element_text(size=20),
           plot.title=element_text(size=30,hjust=0.5),
           panel.grid=element_blank()
   )
   保存图像：
   > ggsave("freq_plot.pdf",plot=freq_plot,width=8,height=8)
   
   
   图2. SNP等位基因频率变异与环境变量间相关联的频率分布直方图
   
   以log₁₀[Bayes factor (BF)] > 1.5为标准，可统计与特定环境变量强相关的SNP数目，这些SNPs通常被认为与群体的本地适应相关。在实践中，环境关联分析可与F_ST离群值分析 (详见本手册另一篇文章 (金铃和郭宝成，2021)) 等结合进行。
   

致谢

我们的研究受到国家自然科学基金委员会优秀青年科学基金项目 (32022009)、中国科学院率先行动百人计划项目和第二次青藏高原综合科学考察研究项目 (Grant No. 2019QZKK0501) 的资助支持。

竞争性利益声明

作者声明没有利益冲突。

参考文献

1. 金铃，郭宝成. (2021). 群体基因组水平的适应性位点检测. Bio-101: e1010627.
2. Ahrens, C. W., Rymer, P. D., Stow, A., Bragg, J., Dillon, S., Umbers, K. D. L. and Dudaniec, R. Y. (2018). The search for loci under selection: trends, biases and progress. Mol Ecol 27(6): 1342-1356.
3. Coop, G., Witonsky, D., Di Rienzo, A. and Pritchard, J. K. (2010). Using environmental correlations to identify loci underlying local adaptation. Genetics 185(4): 1411-1423.
4. Fraik, A. K., Margres, M. J., Epstein, B., Barbosa, S., Jones, M., Hendricks, S., Schönfeld, B., Stahlke, A. R., Veillet, A., Hamede, R., McCallum, H., Lopez-Contreras, E., Kallinen, S. J., Hohenlohe, P. A., Kelley, J. L. and Storfer, A. (2020). Disease swamps molecular signatures of genetic-environmental associations to abiotic factors in Tasmanian devil (Sarcophilus harrisii) populations. Evolution 74(7): 1392-1408.
5. Gunther, T. and Coop, G. (2013). Robust identification of local adaptation from allele frequencies. Genetics 195(1): 205-220.
6. Guo, B., DeFaveri, J., Sotelo, G., Nair, A. and Merilä, J. (2015). Population genomic evidence for adaptive differentiation in Baltic Sea three-spined sticklebacks. BMC Biol 13(1): 19.
7. Guo, B., Li, Z. and Merilä, J. (2016a). Population genomic evidence for adaptive differentiation in the Baltic Sea herring. Mol Ecol 25(12): 2833-2852.
8. Guo, B., Lu, D., Liao, W. B. and Merilä, J. (2016b). Genomewide scan for adaptive differentiation along altitudinal gradient in the Andrew's toad Bufo andrewsi. Mol Ecol 25(16): 3884-3900.
9. Jeffreys, H. (1961). Theory of probability. Oxford University Press. London. ISBN: 0-19-850368-7.
10. Liberles, D. A., Chang, B., Geiler-Samerotte, K., Goldman, A., Hey, J., Kacar, B., Meyer, M., Murphy, W., Posada, D. and Storfer, A. (2020). Emerging frontiers in the study of molecular evolution. J Mol Evol 88(3): 211-226.
11. Magalhaes, I. S., Whiting, J. R., D'Agostino, D., Hohenlohe, P. A., Mahmud, M., Bell, M. A., Skulason, S. and MacColl, A. D. C. (2021). Intercontinental genomic parallelism in multiple three-spined stickleback adaptive radiations. Nat Ecol Evol 5(2): 251-261.
12. Pardo-Diaz, C., Salazar, C. and Jiggins, C. D. (2015). Towards the identification of the loci of adaptive evolution. Methods Ecol Evol 6(4): 445-464.
13. Rellstab, C., Gugerli, F., Eckert, A. J., Hancock, A. M. and Holderegger, R. (2015). A practical guide to environmental association analysis in landscape genomics. Mol Ecol 24(17): 4348-4370.
    
    

附表1. 可进行环境关联分析的方法和软件概览 其中一些方法可用其他软件或R包实现。(译自Rellstab et al., 2015 Table 1)