Evolutionary drivers of sexual signal variation in Amazon slender anoles

Phenotypic variation among populations, as seen in the signaling traits of many species, provides an opportunity to test whether similar factors generate repeated phenotypic patterns in different parts of a species’ range. We investigated whether genetic divergence, abiotic gradients, and sympatry with closely related species explain variation in the dewlap colors of Amazon Slender Anoles, Anolis fuscoauratus. To this aim, we characterized dewlap diversity in the field with respect to population genetic structure and evolutionary relationships, assessed whether dewlap phenotypes are associated with climate or landscape variables, and tested for non-random associations in the distributions of A. fuscoauratus phenotypes and sympatric Anolis species. We found that dewlap colors vary among but not within sites in A. fuscoauratus. Regional genetic clusters included multiple phenotypes, while populations with similar dewlaps were often distantly related. Phenotypes did not segregate in environmental space, providing no support for optimized signal transmission at a local scale. Instead, we found a negative association between certain phenotypes and sympatric Anolis species with similar dewlap color attributes, suggesting that interactions with closely related species promoted dewlap divergence among A. fuscoauratus populations. Amazon Slender Anoles emerge as a promising system to address questions about parallel trait evolution and the contribution of signaling traits to speciation.

Methods
Brief description of the methods for data collection:

Dewlap color variation data

To characterize geographic dewlap color variation in Anolis fuscoauratus, we used data from our comprehensive herpetofaunal inventories in Amazonia and the Atlantic Forest over the last two decades. To this purpose, we sampled individuals by hand or pitfall traps. We only included in our environmental and species co-occurrence analyses 32 sites for which dewlap color information was available (pictures or field notes) from the 63 sites that were included in genetic analyses.

We also obtained data on the presence of all other anole species sympatric with A. fuscoauratus at a given site based on our field inventory data. To reduce the chance of undetected species, we only included data from surveys that lasted a minimum of one week and involved at least three herpetologists searching for animals both night and day. The final dataset included occurrence data for 11 other anole species at the 32 sites for which A. fuscoauratus dewlap coloration data were available: Anolis auratus, A. chrysolepis, A. dissimilis, A. nasofrontalis, A. ortonii, A. planiceps, A. punctatus, A. scypheus, A. tandai, A. trachyderma, and A. transversalis.

Genetic data

A double-digest restriction site associated DNA library (ddRAD) was generated at the University of Wisconsin Biotechnology Center. Briefly, DNA extractions were digested with the restriction enzymes PstI and MspI, and the resulting fragments were tagged with individual barcodes, PCR-amplified, multiplexed, and sequenced in a single lane on an Illumina HiSeq 2500 platform. The number of paired-end reads ranged from ~1.15 to 8.85 million per individual, with a read length of 100 base pairs. De-multiplexed raw sequence data were deposited in the Sequence Read Archive (BioProject PRJNA492310; BioSample accessions SAMN18340748-18340924).

We used Ipyrad v. 0.7.30 to de-multiplex and assign reads to individuals based on sequence barcodes (allowing no mismatches from individual barcodes), perform de novo read assembly (minimum clustering similarity threshold = 0.95), align reads into loci, and call single nucleotide polymorphisms (SNPs). A minimum Phred quality score (= 33), sequence coverage (= 6x), read length (= 35 bp), and maximum proportion of heterozygous sites per locus (= 0.5) were enforced, while ensuring that variable sites had no more than two alleles (i.e., a diploid genome). Moreover, for inclusion in the final datasets, we ensured that each locus was present in at least 70% of the sampled individuals.

To estimate population genetic structure and admixture in Anolis fuscoauratus, we generated in Ipyrad a final dataset composed of 118,434 SNPs at 16,368 loci (including no outgroups). A single SNP was then extracted from each locus to minimize sampling of linked SNPs. We used VCFtools v. 0.1.16 to filter out SNPs whose minor allele frequency (MAF) was lower than 0.05.

To estimate phylogenetic relationships, we generated in Ipyrad a second dataset composed of 135,952 SNPs at 17,302 RAD loci (now including outgroup taxa and linked SNPs), ensuring that each locus was present in at least 70% of the sampled individuals.

Environmental data

We used 17 variables in environmental analyses: cover of evergreen broadleaf trees, deciduous broadleaf trees, shrubs, herbaceous vegetation, and regularly flooded vegetation, annual cloud cover, elevation, slope, terrain roughness, and terrain ruggedness, all obtained from the EarthEnv database. As climatic variables, we used annual mean temperature, maximum temperature of the warmest month, mean temperature of the warmest quarter, annual precipitation, precipitation of the wettest month, and precipitation of the wettest quarter, obtained from the Chelsa database, as well as the climatic moisture index, a metric of relative wetness, obtained from the ENVIREM database. Values were extracted for each environmental variable from the 32 sites for which A. fuscoauratus dewlap color information was available in QGIS v. 3.4.5.

Please refer to the manuscript's Material and Methods for details on the analyses.

Please send questions to ivanprates [at] gmail [dot] com

诸多物种种群间的表型变异（phenotypic variation）——如诸多物种的信号特征所呈现的那样——为检验相似因子是否在物种分布的不同区域反复催生了一致的表型模式提供了研究契机。本研究以亚马逊细安乐蜥（Anolis fuscoauratus）的喉扇（dewlap）色彩变异为研究对象，探究遗传分化、非生物梯度以及与近缘物种同域分布（sympatry）能否解释该变异。

为此，我们基于野外调查对喉扇多样性进行了系统表征，结合种群遗传结构与演化关系展开分析，评估了喉扇表型（phenotype）是否与气候或景观变量存在关联，并检验了亚马逊细安乐蜥的表型分布与同域安乐蜥物种是否存在非随机关联。

研究结果显示，亚马逊细安乐蜥的喉扇色彩仅存在种群间差异，而同一生境内无显著变异。区域遗传聚类涵盖多种表型，而喉扇特征相似的种群往往亲缘关系较远。表型并未在环境空间中出现分化，未支持局部尺度下的最优信号传递假说。与之相反，我们发现部分表型与拥有相似喉扇色彩特征的同域安乐蜥物种呈负相关，提示与近缘物种的相互作用推动了亚马逊细安乐蜥种群间的喉扇分化。本研究使亚马逊细安乐蜥成为探究平行性状演化以及信号特征对物种形成贡献相关问题的极具潜力的研究系统。

## 研究方法

### 数据采集方法简要说明
#### 喉扇色彩变异数据
为表征亚马逊细安乐蜥的地理喉扇色彩变异，我们使用了过去二十年间在亚马逊流域与大西洋森林开展的全面两栖爬行动物调查数据。我们通过手捕或陷阱诱捕的方式采集个体。本研究的环境与物种共存分析仅纳入了63个遗传分析位点中的32个具备喉扇色彩信息（照片或野外记录）的采样点。

我们还基于野外调查数据，获取了各采样点处与亚马逊细安乐蜥同域分布的所有其他安乐蜥物种的存在数据。为降低未检测到物种的概率，我们仅纳入了持续至少一周、且至少有3名两栖爬行动物学家昼夜开展搜索的调查数据。最终数据集纳入了具备亚马逊细安乐蜥喉扇色彩数据的32个采样点处的11种其他安乐蜥物种的出现记录：金项安乐蜥（Anolis auratus）、金鳞安乐蜥（Anolis chrysolepis）、异色安乐蜥（Anolis dissimilis）、鼻额安乐蜥（Anolis nasofrontalis）、奥顿安乐蜥（Anolis ortonii）、平头安乐蜥（Anolis planiceps）、斑点安乐蜥（Anolis punctatus）、舟形安乐蜥（Anolis scypheus）、坦达伊安乐蜥（Anolis tandai）、粗皮安乐蜥（Anolis trachyderma）以及横带安乐蜥（Anolis transversalis）。

#### 遗传数据
我们在威斯康星大学生物技术中心构建了双酶切限制性位点关联DNA文库（double-digest restriction site associated DNA library，ddRAD）。简要而言，我们使用限制性内切酶PstI与MspI对DNA提取物进行酶切，所得片段标记有个体特异性条形码（barcode），经PCR扩增、多重化后，在Illumina HiSeq 2500平台的单个泳道中完成测序。每个个体的双端读段（paired-end reads）数量约为115万至885万，读段长度为100碱基对。去多路复用的原始序列数据已提交至序列读取档案库（Sequence Read Archive，BioProject编号PRJNA492310；生物样本登录号SAMN18340748-18340924）。

我们使用Ipyrad v. 0.7.30对序列进行去多路复用，并依据序列条形码将读段分配至对应个体（允许个体条形码无错配），进行从头读段组装（最小聚类相似性阈值=0.95），将读段比对至基因座并调用单核苷酸多态性（single nucleotide polymorphisms，SNPs）。分析过程中设置了以下过滤参数：最小Phred质量分数=33、最小序列覆盖度=6倍、最小读段长度=35碱基对、每个基因座的杂合位点最大比例=0.5，同时确保变异位点仅存在两个等位基因（即二倍体基因组）。此外，为纳入最终数据集，我们要求每个基因座在至少70%的采样个体中存在。

为估计亚马逊细安乐蜥的种群遗传结构与基因混合（admixture），我们在Ipyrad中生成了包含16368个基因座上的118434个SNPs的最终数据集（不包含外类群）。随后从每个基因座中提取单个SNP，以降低连锁SNPs的采样偏差。我们使用VCFtools v. 0.1.16过滤掉次要等位基因频率（minor allele frequency，MAF）低于0.05的SNPs。

为估计系统发育关系，我们在Ipyrad中生成了第二个数据集，包含17302个RAD基因座上的135952个SNPs（此时包含外类群类群与连锁SNPs），并要求每个基因座在至少70%的采样个体中存在。

#### 环境数据
我们在环境分析中使用了17个变量：常绿阔叶树覆盖率、落叶阔叶树覆盖率、灌木覆盖率、草本植被覆盖率、常淹植被覆盖率、年云量、海拔、坡度、地形粗糙度与地形崎岖度，所有数据均取自EarthEnv数据库。气候变量方面，我们使用了来自Chelsa数据库的年平均气温、最热月最高气温、最热季度平均气温、年降水量、最湿月降水量与最湿季度降水量，以及来自ENVIREM数据库的气候湿度指数（衡量相对湿润程度的指标）。我们使用QGIS v. 3.4.5从32个具备亚马逊细安乐蜥喉扇色彩数据的采样点中提取了各环境变量的数值。

有关分析的详细细节，请参阅论文的“材料与方法”部分。

如有疑问，请致信ivanprates [at] gmail [dot] com