DOI: 10.5061/dryad.r2280gbnh
title: The complex interplay between evolutionary flexibility and diversification in a family of spiders

description: Understanding the mechanisms of how morphological traits drive speciation
 and contribute to species richness is pivotal in evolutionary biology. In
 this context, the evolutionary flexibility of morphological traits may
 play a significant role. Using the diverse daddy long-legs spiders,
 Pholcidae, which currently includes almost 2,000 described species, we
 explore the interplay between speciation rate, trait evolution rate,
 microhabitat shift rate, species richness, interspecific variability of
 body size, leg length, relative leg length, and leg proportions. We apply
 a combination of large-scale genomic and taxonomic sampling, and
 phylogenetic and comparative analyses to assess the dynamics of
 diversification and evolutionary flexibility (measured as either the
 standard variance or disparity of traits), as well as their interactions.
 We found that increased evolutionary flexibility is accompanied by
 accelerated rates in speciation and trait evolution, and by higher species
 richness. We also observed near-isometry of leg length and body size in
 the species-rich lineages, suggesting stabilizing selection, while
 positive allometry in the species-poor lineages indicates directional
 selection. Additionally, we found a positive correlation between trait
 evolution rates and microhabitat shifts. We argue that environmental
 heterogeneity and frequent microhabitat shifts may contribute to the
 origin and maintenance of evolutionary flexibility, which in turn
 influences the organisms’ ability to exploit new resources and habitats.
 On the other hand, our study suggests that a lack of evolutionary
 flexibility, e.g., due to dwarfism or gigantism, could be an evolutionary
 dead end. Our findings also support the idea that taxa do not always
 follow Cope’s rule, which postulates that body size increases with lineage
 age. This study enhances our understanding of the mechanisms of
 diversification, demonstrating that the evolutionary flexibility of
 morphological traits may vary among closely related taxa and is likely to
 have contributed to the uneven distribution of biodiversity across the
 tree of life.
description: # The complex interplay between evolutionary flexibility and
 diversification in a family of spiders URL:
 [https://doi.org/10.5061/dryad.r2280gbnh](https://doi.org/10.5061/dryad.r2280gbnh) **Authors**: Guanliang Meng^1^, Lars Podsiadlowski^1^, Dimitar Dimitrov^2^, and Bernhard A. Huber^1^ 1 Zoological Research Museum Alexander Koenig, LIB, Bonn, Germany 2 University Museum of Bergen, University of Bergen, Bergen, Norway Corresponding author: Guanliang Meng ([G.Meng@leibniz-lib.de](mailto:G.Meng@leibniz-lib.de)); Dimitar Dimitrov ([dimitard.gwu@gmail.com](mailto:dimitard.gwu@gmail.com)) **ORCiD**: Guanliang Meng: [https://orcid.org/0000-0002-6488-1527](https://orcid.org/0000-0002-6488-1527) Lars Podsiadlowski: [https://orcid.org/0000-0001-7786-8930](https://orcid.org/0000-0001-7786-8930) Dimitar Dimitrov: [https://orcid.org/0000-0001-5830-5702](https://orcid.org/0000-0001-5830-5702) Bernhard A. Huber: [https://orcid.org/0000-0002-7566-5424](https://orcid.org/0000-0002-7566-5424) ## 1. SUPPLEMENTARY MATERIAL 1\. **Supplementary Text and Figures**: `Supplementary_Text_and_Figures.pdf` 2\. **Supplementary Tables**: `Supplementary_Tables.zip` 3\. **Data (e.g., sequences, alignments) and scripts for target capture sequencing and comparative analysis**: `data_and_scripts.zip` 4\. **UCE+Sanger-data-based trees**: `Supplementary_Data_2.zip` 5\. **PDF version of the README file**: `README.pdf`. ### 2. Further details on the file `data_and_scripts.zip` This file contains the data, scripts, and commands used by the present study. Folder structure: ``` . ├── 1_script │   └── FASconCAT-G_v1.05.pl ├── 2_phyluce │   ├── 1_raw_data │   ├── 2_clean_data │   ├── 3_assembly │   ├── 4_detect_cross_contamination │   ├── 5_match_contigs_to_probes │   └── 6_extract_uce_loci ├── 3_UCE_phylogeny │   ├── 1_trimmed_msa │   ├── 2_tree_construction │   └── 3_dating ├── 4_long-name-short-name-rela │   ├── June-30-2024.six-genes-samples.short2long.tsv │   └── June-30-2024.UCE-samples.short2long.tsv ├── 5_six_gene_phylogeny │   ├── 1_trimmed_msa │   ├── 2_tree_construction │   └── 3_dating ├── 6_comparative_analysis │   ├── 0.dated_trees │   ├── 1.diversification.BAMM │   ├── 2.trait_evolution.BAMM.subfam-level-sampling-fraction │   ├── 3.model_fitting │   ├── 4.phylogenetic_signal │   ├── 5.ancestral_state │   ├── 6.relationships │   ├── 7.disparity.28-05-2024 │   ├── 8.linage_through_time │   ├── readme.txt │   └── Trait data for Meng et al 2025.xlsx ├── 7_UCE_and_six-gene_phylogeny │   ├── 1_unpartitioned │   ├── 2_partitioned │   ├── April-9-2024.code.subfam.tsv │   ├── combined_UCE+Sanger_trees.pdf │   └── combined.code.subfam.tsv └── readme.txt ``` ## 2.1 Folder `1_script/` This folder contains the Perl script `FASconCAT-G_v1.05.pl`, which was from the publication: Kück, P., Longo, G.C. FASconCAT-G: extensive functions for multiple sequence alignment preparations concerning phylogenetic studies. *Front Zool* **11**, 81 (2014). [https://doi.org/10.1186/s12983-014-0081-x](https://doi.org/10.1186/s12983-014-0081-x) ## 2.2 Folder `2_phyluce/` This folder contains the commands for generating ultra-conserved elements (UCEs) from raw sequencing data using the Phyluce pipeline (Faircloth 2016). This folder has the following subfolders: ``` ├──1_raw_data ├──2_clean_data ├──3_assembly ├──4_detect_cross_contamination ├──5_match_contigs_to_probes └──6_extract_uce_loci ``` ### 2.2.1 Subfolder `1_raw_data/` Please download the raw data from the NCBI SRA database (BioProject No: PRJNA851037), see Supplementary Table S1 for more details. ### 2.2.2 Subfolder `2_clean_data/` **Task**: Raw data filtering. This folder contains the command for running fastp (v0.23.1) (Chen et al. 2018), which is to filter raw data from the folder `1_raw_data` and generate clean data. ​ Briefly, reads would be filtered out if: (1) one read contained more than ten Ns; (2) more than 40% of bases were unqualified; (3) reads were duplicates; (4) reads were shorter than 40 bp. ### 2.2.3 Subfolder `3_assembly/` **Task**: UCE contig assembly. **Description**: This folder contains the command for running the executable `phyluce_assembly_assemblo_spades` from the Phyluce pipeline (Faircloth 2016). The result is a contig file for each sample. Default parameters were used except for `--cores 16 --memory 100` (i.e., 16 CPUs and a maximum of 100 GB RAM would be applied) The file `assembly.conf` contains the "code" (a short name for each sample) and the path to the folder containing fastq 1 and fastq 2 files of each sample. Users should adjust this accordingly. When the assembly for all samples is finished, we may copy all assembled contig files into the same directory (e.g., `spades-assemblies/contigs`) ### 2.2.4 Subfolder `4_detect_cross_contamination/` **Task**: Filter out potential cross-sample contaminants using Blast. To this end, we firstly make each sequence ID unique by adding the sample "code" to the sequence ID (see the folder `1_add_sampleID_to_seq-title/`), before merging all sequences of all samples into a single file (`2_all-blast-all/all_assemblies.fa`). Next, the `2_all-blast-all/all_assemblies.fa` file was blasted against itself using the blastn command with default parameters. Blast results with "query" and "subject" coming from the same sample are then excluded. Blast results with alignment similarity &gt;= 99% and alignment length &gt;= 100 bp are considered as potential contamination. However, they could also be highly conserved elements such as rDNA genes. Finally, we used the following criteria to keep or remove potential contaminants (see the folder `3_filter_based_on_abundance/`): if the abundance difference (output by SPAdes) between the query and subject sequences was less than 10X, both sequences were marked as 'removed'. If both sequences had abundances ≥ 20X and their abundance difference exceeded 10X, the sequence with lower abundance was marked as 'removed'. Any sequences marked as 'removed' were excluded from the assembly. All the clean sequences are presented in the folder `4_clean_assemblies/`. ### 2.2.5 Subfolder `5_match_contigs_to_probes/` **Task**: Match the assembled contigs to the probe set (see the file`SpiderUCE_probe_set_men13099-sup-0002-Supinfo.fasta`) using the command `phyluce_assembly_match_contigs_to_probes` from the Phyluce pipeline. The probe set is provided by the following publication: Kulkarni, S., Wood, H., Lloyd, M., and Hormiga, G. (2020), Spider-specific probe set for ultraconserved elements offers new perspectives on the evolutionary history of spiders (Arachnida, Araneae), Mol. Ecol. Resour., 20 (1), 185-203. ### 2.2.6 Subfolder `6_extract_uce_loci/` **Task**: Extact homologous UCE among species using the command `phyluce_assembly_get_fastas_from_match_counts` from the Phyluce pipeline. Matrices 1 to 10 with various gene and taxon coverages are created, with the help of PhyloMAd, BMGE, and RY-coding. They are placed in the folder `../3_UCE_phylogeny/1_trimmed_msa/`. ## 2.3 Folder `3_UCE_phylogeny/` **Task**: Phylogenetic inference using UCE datasets. This folder has the following subfolders: ``` ├──1_trimmed_msa ├──2_tree_construction └──3_dating ``` ### 2.3.1 Subfolder `1_trimmed_msa/` This folder contains the sequences of Matrices 1 to 10. Each matrix has its own folder. ``` ├── matrix_1 ├── matrix_2 ├── matrix_3 ├── matrix_4 ├── matrix_5 ├── matrix_6 ├── matrix_7 ├── matrix_8 ├── matrix_9 └── matrix_10 ``` To create a concatenated alignment for each matrix, users should enter each directory and run the following command, for instance, ```bash # concatenate the alignments to Phylip format cd 3_UCE_phylogeny/1_trimmed_msa/matrix_3 perl ../../../1_script/FASconCAT-G_v1.05.pl -s -p -p -n -l ``` The files named `FcC_supermatrix.phy` and `FcC_supermatrix_partition.txt` will then be created, which can be used for subsequent tree construction. ### 2.3.2 Subfolder `2_tree_construction/` This folder contains tree construction commands and the **resulting trees** for each matrix. Each matrix has its own folder: ``` . ├── matrix_1 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-1.gblocks.IQ-TREE.partitioned.tre ├── matrix_10 │   └── 1_unPartitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-10.IQ-TREE.unpartitioned.tre ├── matrix_2 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-2.gblocks.IQ-TREE.partitioned.tre ├── matrix_3 │   ├── 1_Partitioned_ML.IQ-TREE │   │   ├── iqtree.sh │   │   └── Matrix-3.TrimAl.IQ-TREE.partitioned.clean-support.tre │   ├── 2_Partitioned_ML.RAxML-NG │   │   ├── Matrix-3.TrimAl.RAxML-NG.partitioned.tre │   │   └── raxml.sh │   ├── 3_PhyloBayes-MPI │   │   ├── check │   │   │   └── record.sh │   │   ├── Matrix-3.TrimAl.phylogbayes-MPI.unpartitioned.tre │   │   ├── run1 │   │   │   ├── b.sh │   │   │   └── sub.sh │   │   └── run2 │   │   ├── b.sh │   │   └── sub.sh │   ├── 4_ASTRAL │   │   ├── 1.locus_tree │   │   │   └── batch_tree.sh │   │   └── 2_species_tree │   │   ├── astral.sh │   │   └── Matrix-3.TrimAl.Astral.UFBS30.clean-support.tre │   ├── 5_SVDquartets │   │   ├── 1_individual_msa │   │   ├── Matrix-3.TrimAl.SVDquartets.unpartitioned.tre │   │   ├── record.sh │   │   └── run1 │   │   ├── cmd1.txt │   │   └── run1.sh │   └── 6_Constrained_partitioned │   ├── monophyletic_other_Arteminae_plus_Artema │   │   ├── Arteminae.tre │   │   └── multi-partition.iqtree.sh │   └── Ninetinae_sister_to_all_other_pholcids │   ├── multi-partition.iqtree.sh │   └── Ninetaine_sister_to_all_others.tre ├── matrix_4 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-4.TrimAl.IQ-TREE.partitioned.tre ├── matrix_5 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-5.TrimAl.IQ-TREE.partitioned.tre ├── matrix_6 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-6.TrimAl.IQ-TREE.partitioned.tre ├── matrix_7 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-7.IQ-TREE.partitioned.tre ├── matrix_8 │   └── 1_Partitioned_ML.IQ-TREE │   ├── iqtree.sh │   └── Matrix-8.IQ-TREE.partitioned.tre └── matrix_9 └── 1_unPartitioned_ML.IQ-TREE ├── iqtree.sh └── Matrix-9.IQ-TREE.unpartitioned.tre ``` #### Example I will explain the tree construction steps using `matrix_3/` as an example. The folder `matrix_3/` has the following content: ``` ├── 1_Partitioned_ML.IQ-TREE ├── 2_Partitioned_ML.RAxML-NG ├── 3_PhyloBayes-MPI ├── 4_ASTRAL ├── 5_SVDquartets └── 6_Constrained_partitioned ``` **IQ-TREE** In the folder `1_Partitioned_ML.IQ-TREE/`, we perform partitioned analysis using the maximum likelihood (ML) method IQ-TREE: ```bash # Partitioned ML analysis with IQ-TREE aln=../../../1_trimmed_msa/matrix_3/FcC_supermatrix.phy partition=../../../1_trimmed_msa/matrix_3/FcC_supermatrix_partition.txt /home/gmeng/soft/iqtree-2.2.0-Linux/bin/iqtree2 \ -s $aln \ -p $partition \ --prefix m3_partitioned \ -m MFP+MERGE \ --sampling GENESITE \ -pers 0.2 -nstop 500 -B 2000 -bnni --alrt 2000 -T 24 --runs 10 &gt;i.log 2&gt;i.err ``` Best-fit models of molecular evolution were determined by ModelFinder (Kalyaanamoorthy et al. 2017) in IQ-TREE automatically (-m MFP+MERGE). To overcome local optima during heuristics, we performed ten independent runs on each dataset, with a smaller perturbation strength (-pers 0.2) and a larger number of stop iterations (-nstop 500) (Nguyen et al. 2015; Liu et al. 2024). Branch supports were evaluated with 2000 ultrafast bootstrap (UFBoot) (Minh et al. 2013) considering potential model violations (-B 2000 -bnni --sampling GENESITE). The SH-aLRT branch test (Guindon et al. 2010) was performed using 2,000 bootstrap replicates. **RAxML-NG** In the folder `2_Partitioned_ML.RAxML-NG/`, we perform ML analysis using RAxML-NG. ```bash aln=../../../1_trimmed_msa/matrix_3/FcC_supermatrix.phy partition=../../../1_trimmed_msa/matrix_3/FcC_supermatrix_partition.txt sed 's#DNA#GTR+R4+I+FO#g' $partition &gt; partition.txt # RAxML-NG v. 1.1 released on 29.11.2021 by The Exelixis Lab. /home/gmeng/.conda/envs/mybase/bin/raxml-ng \ --all \ --threads 8 \ --msa $aln \ --model partition.txt \ --brlen scaled \ --bs-trees autoMRE \ --tree pars{25},rand{25} \ --seed 2 \ --prefix T1 ``` During each RAxML-NG run, the GTR+R4+I+FO model, and 25 parsimony starting trees + 25 random starting trees were used, and the number of bootstrap replicates was automatically determined (--bs-trees autoMRE) (Pattengale et al. 2009). **PhyloBayes-MPI** In the folder `3_PhyloBayes-MPI/`, to model both rate-heterogeneity across sites and pattern-heterogeneity across sites (i.e., across-site compositional heterogeneity), we applied the CAT-GTR+G4 model (Lartillot and Philippe 2004) in PhyloBayes-MPI. It is less sensitive to systematic errors than classical empirical models (Lartillot et al. 2013). ```bash module load phylobayes_mpi/1.8 export OMPI_MCA_opal_cuda_support=false aln=../../../../1_trimmed_msa/matrix_3/FcC_supermatrix.phy mpirun -np $NSLOTS pb_mpi -d $aln -cat -gtr -dgam 4 T1 ``` Note: At least two runs (T1, T2) should be performed, and MCMC chain convergence must be examined: ```bash module load phylobayes_mpi/1.8 export OMPI_MCA_opal_cuda_support=false tracecomp -x 20000 T1.trace T2.trace &gt; tracecomp.log 2&gt;&1 bpcomp -x 20000 T1.treelist T2.treelist &gt; bpcomp.log ``` Please check the manual of PhyloBayes-MPI to learn how to determine the convergence of MCMC chains. **ASTRAL** Unlike the previous concatenation-based method, to account for incomplete lineage sorting (ILS) and mitigate the impact of missing data, we perform two-step multi-species coalescence (MSC) tree inference using ASTRAL in the folder `4_ASTRAL/`. During the process, gene trees are first inferred using IQ-TREE: ```bash cd 3_UCE_phylogeny/2_tree_construction/matrix_3/4_ASTRAL/1.locus_tree ls ../../../../1_trimmed_msa/matrix_3/*fasta | while read f ; do a=$(basename $f .fasta) ; mkdir -p $a ; echo "/home/gmeng/soft/iqtree-2.2.0-Linux/bin/iqtree2 -s ../$f --prefix $a -pers 0.2 -nstop 500 -B 2000 -bnni -T 2 &gt;i.log 2&gt;i.err" &gt;$a/loci.iqtree.sh ; cd $a ; sh loci.iqtree.sh ; cd ../ ; done ``` The internal nodes with low ultra-fast bootstrap supports (&lt;= 30) of gene trees are collapsed, and after which, a species tree is summarized: ```bash cd 3_UCE_phylogeny/2_tree_construction/matrix_3/4_ASTRAL/2_species_tree # cat all genes cat ../1.locus_tree/*/*.treefile &gt;gene_trees.txt nw_ed gene_trees.txt 'i & b&lt;=30' o &gt; gene_trees.UFBS30.txt java_lib=/home/gmeng/soft/ASTRAL-single-thread-v5.7.8/Astral/ jar=/home/gmeng/soft/ASTRAL-single-thread-v5.7.8/Astral/astral.5.7.8.jar java \ -Djava.library.path=$java_lib \ -jar $jar \ -i gene_trees.UFBS30.txt \ -o species.UFBS30.tre \ 2&gt;species.UFBS30.tre.err ``` **SVDquartets** The ASTRAL method is sensitive to gene tree error, for instance, due to the few informative sites in the short alignments. Therefore, we also tested for the SVDquartets, which infers a species tree based on site patterns based on MSC. The parameters file ` cmd1.txt`: ```bash log file=q1.log; svdq evalQuartets=all seed=135454568 treeInf=QFM bootstrap=multilocus loci=combined nreps=200 nthreads=24 treemodel=mscoalescent treefile=concat.SVDquartets.bs.tre replace=yes ambigs=missing plus2=yes ; saveTrees file=SVDquartets.run3.species.tre supportValues=nodeLabels ; ``` Run the analysis: ```bash /home/gmeng/soft/bin/paup4a168_centos64 ../concat_msa_nexus/concat_msa_nexus.loci-renamed.nexus &lt; cmd1.txt ``` You may test different seeds and check if the topologies of the species trees are consistent. **Constrained tree search and topology test** We tested two hypotheses: (1) the traditional Arteminae (i.e., ‘other Arteminae’ + *Artema*) is a monophyletic group. (2) Ninetinae is sister to all other pholcids. ``` . ├── monophyletic_other_Arteminae_plus_Artema │   ├── Arteminae.tre │   └── multi-partition.iqtree.sh └── Ninetinae_sister_to_all_other_pholcids ├── multi-partition.iqtree.sh └── Ninetaine_sister_to_all_others.tre ``` ### 2.3.3 Subfolder `3_UCE_phylogeny/3_dating/` **Task**: Divergence time calibration using MCMCTREE. We followed Magalhaes et al. (2020) for the justifications of the fossil assignments; see Supplementary Table S3 for details. The oldest representative fossil of each clade was used to calibrate its divergence from the sister group (stem-dating), following Benton and Donoghue (2006) and Magalhaes et al. (2020). The soft-bound strategy was employed. We set the age bounds of the root (Synspermiata + Hypochilidae) as follows: the minimum bound (164 million years ago, Ma) reflects the placement of the genus *Eoplectreurys* as a stem-Synspermiata, based on the oldest known fossil species, *E. gertschi* Selden & Huang (Gao and Ren 2006; Magalhaes et al. 2020); the conservative maximum bound (400 Ma) was chosen because no Araneae stem fossils older than 383 Ma are known (Magalhaes et al. 2020). Two independent runs are performed, and the convergence of MCMC chains is examined. For learning how to perform divergence time calibration using MCMCTREE, please refer to: * Dos Reis M, Yang Z. Bayesian Molecular Clock Dating Using Genome-Scale Datasets. Methods Mol Biol. 2019;1910:309-330. doi: 10.1007/978-1-4939-9074-0_10. PMID: 31278669. * [https://github.com/mariodosreis/divtime](https://github.com/mariodosreis/divtime) ## 2.4 Folder `4_long-name-short-nam/` **Description**: These two files below indicate the relationships between the sample "code" and long sample names. ``` . ├── June-30-2024.six-genes-samples.short2long.tsv └── June-30-2024.UCE-samples.short2long.tsv ``` ## 2.5 Folder `5_six_gene_phylogeny/` **Task**: Phylogenetic tree inference using the six-gene matrices + divergence time calibration. This folder has the following content: ``` . ├── 1_trimmed_msa ├── 2_tree_construction └── 3_dating ``` ### 2.5.1 Subfolder `1_trimmed_msa`/ This folder contains alignments derived from 10 different trimming strategies: ``` . ├── trimmed_with_ClipKIT_gappy ├── trimmed_with_ClipKIT_kpi ├── trimmed_with_ClipKIT_kpi-gappy ├── trimmed_with_ClipKIT_kpi-smart-gap ├── trimmed_with_ClipKIT_kpic ├── trimmed_with_ClipKIT_kpic-gappy ├── trimmed_with_ClipKIT_kpic-smart-gap ├── trimmed_with_ClipKIT_smart-gap ├── trimmed_with_Gblocks └── trimmed_with_TrimAl ``` ### 2.5.2 Subfolder `2_tree_construction/` **Task**: To perform phylogenetic tree inference using the six-gene matrices. As a constraint (the file `clean.constraint.tre`), we used the tree from the IQ-TREE partitioned ML analysis of matrix 3. This folder contains the tree construction commands and the **resulting trees**: ``` . ├── trimmed_with_ClipKIT_gappy │   ├── clean.constraint.tre │   ├── clean.input.ultrametric.tree │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_kpi │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.ClipKIT_kpi.partitioned.IQ-TREE.rooted.treePL-dated.tre │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_kpi-gappy │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.ClipKIT_kpi-gappy.partitioned.IQ-TREE.rooted.treePL-dated.tre │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_kpi-smart-gap │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.ClipKIT_kpi-smart-gap.partitioned.IQ-TREE.rooted.treePL-dated.tre │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_kpic │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.ClipKIT_kpic.partitioned.IQ-TREE.rooted.treePL-dated.tre │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_kpic-gappy │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.ClipKIT_kpic-gappy.partitioned.IQ-TREE.rooted.treePL-dated.tre │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_kpic-smart-gap │   ├── clean.constraint.tre │   ├── clean.input.ultrametric.tree │   └── multi-partition.iqtree.sh ├── trimmed_with_ClipKIT_smart-gap │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.ClipKIT_smart-gap.partitioned.IQ-TREE.rooted.treePL-dated.no-outgroup-and-duplicate.no-outgroup-and-duplicate.binary.ultrametric.tre │   └── multi-partition.iqtree.sh ├── trimmed_with_Gblocks │   ├── clean.constraint.tre │   ├── July-9-2023-v1.six-gene.Gblocks.partitioned.IQ-TREE.rooted.treePL-dated.tre │   └── multi-partition.iqtree.sh └── trimmed_with_TrimAl ├── clean.constraint.tre ├── July-9-2023-v1.six-gene.TrimAl.partitioned.IQ-TREE.rooted.treePL-dated.tre └── multi-partition.iqtree.sh ``` ### 2.5.3 Subfolder `3_dating/` **Task**: Calibrating divergence times in the six-gene trees based on the dated UCE-matrix3-IQ-TREE phylogeny. Divergence times of the six-gene trees are estimated using TreePL (v1.0) (Smith and O'Meara 2012) and the R function congruify.phylo() in geiger (v2.0.11) (Pennell et al. 2014), with our dated UCE-matrix 3 tree (see the section Divergence time calibration) used as the reference tree in all analyses. Since our reference tree is a subset of each of the target trees in terms of tip taxa and topology, the concordant nodes between the reference tree and a target tree can be easily determined by the congruify.phylo() function (Eastman et al. 2013). Thus, the divergence times of all nodes in the reference tree serve as calibration points for the target tree, which are passed to TreePL via congruify.phylo() to calibrate the node ages of the target tree (Eastman et al. 2013). ```R setwd('./') library(geiger) # the input six-gene tree target_tree_file= out_tree= target_tree=read.tree(target_tree_file) # the dated UCE-matrix3-IQ-TREE tree ref_tree &lt;- read.tree('../../../reference_tree.mean-time-only.with-CH-codes.tre') res &lt;- congruify.phylo(ref_tree, target_tree, scale='treePL', ncores = 2) write.tree(res$phy, file=out_tree) ``` ## 2.6 Folder `6_comparative_analysis/` **Tasks**: 1. Modeling speciation/extinction dynamics as well as trait evolution 2. Testing for correlations between speciation/extinction and trait evolution We focus on four traits: 1. body size, using male carapace width as a proxy; 2. leg length, using the combined length of male tibia 1 + metatarsus 1 as a proxy, because these can be measured easily and precisely, and they are available from the taxonomic literature; 3. relative leg length, which is the ratio of leg length to body size, and 4. leg proportion, which is the ratio of male metatarsus 1 to tibia 1; it can be used as a proxy for microhabitat preference (i.e., “space”, “leaf”, and “ground”), in which leaf and space living pholcids tend to have larger leg proportion values (Eberle et al. 2018). Note: These trait data (please refer to the file `Trait data for Meng et al 2025.xlsx`) have also been submitted to the **MorphoBank** ([http://morphobank.org/permalink/?P6008](http://morphobank.org/permalink/?P6008)). During the review process of the manuscript, reviewers may access the data by visiting [https://morphobank.org/index.php/LoginReg/form](https://morphobank.org/index.php/LoginReg/form), then enter the project number **P6008** in the email address field and the reviewer login password (**ClipKIT_smart-gap_2025**). Once the manuscript is accepted, the dataset and this repository ([https://doi.org/10.5061/dryad.r2280gbnh](https://doi.org/10.5061/dryad.r2280gbnh)) will be released and made publicly accessible. **Input trees**:\ To explicitly account for the impact of topological uncertainty on our results while ensuring analytical tractability, we selected three representative trees for subsequent comparative analyses: * a primary tree based on the ClipKIT smart-gap trimmed-MSA; * an alternative tree based on the ClipKIT gappy trimmed-MSAs; * another alternative tree based on the kpic-smart-gap trimmed-MSAs. These trees were chosen based on two criteria: * minimum number of species obviously misplaced at the subfamily level; and * monophyly of uncontested genera. This folder contains the following content: ``` . ├── 0.dated_trees ├── 1.diversification.BAMM ├── 2.trait_evolution.BAMM.subfam-level-sampling-fraction ├── 3.model_fitting ├── 4.phylogenetic_signal ├── 5.ancestral_state ├── 6.relationships ├── 7.disparity.28-05-2024 ├── 8.linage_through_time ├── Trait data for Meng et al 2025.xlsx └── readme.txt ``` The following books are a good start to learn comparative analysis: * *Phylogenetic Comparative Methods* (by Luke J. Harmon; [https://lukejharmon.github.io/pcm/pdf/phylogeneticComparativeMethods.pdf](https://lukejharmon.github.io/pcm/pdf/phylogeneticComparativeMethods.pdf) * *Phylogenetic Comparative Methods in R* (by Liam J. Revell and Luke J. Harmon; [http://www.phytools.org/Rbook/](http://www.phytools.org/Rbook/)) ### 2.6.0 Subfolder `0.dated_trees/` **Description**: Tree files and Trait files. For instance, for the ClipKIT_smart-gap six-gene trees: Filename explanations: * `*.body_size_and_leg_length.*`: relative leg length, which is the ratio of leg length to body size. * `*.leg_length.*`: leg length (using the combined length of male tibia 1 + metatarsus 1 as a proxy) * `*.m_carapace_width.*`: body size (using male carapace width as a proxy) * `*.met_to_tibia_ratio.*`: leg proportion (the ratio of male metatarsus 1 to tibia 1; it can be used as a proxy for microhabitat preference (i.e., “space”, “leaf”, and “ground”), in which leaf and space living pholcids tend to have larger leg proportion values (Eberle et al. 2018)) * `*.microhabitat.*`: the raw microhabitat data (only available for a few species) Tree files: ```bash $ ls 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/*.ultrametric.tre 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.body_size_and_leg_length.final.codes.ultrametric.tre 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.leg_length.final.codes.ultrametric.tre 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.m_carapace_width.final.codes.ultrametric.tre 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.met_to_tibia_ratio.final.codes.ultrametric.tre 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.microhabitat.final.codes.ultrametric.tre ``` Trait files: ```bash $ ls 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/*.final.tsv 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.body_size_and_leg_length.final.tsv 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.leg_length.final.tsv 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.m_carapace_width.final.tsv 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.met_to_tibia_ratio.final.tsv 0.dated_trees/2.partitioned/treePL/step5.maybe.non-monophyletic.only_tips_with_traits/ClipKIT_smart-gap/used.code.sp.microhabitat.final.tsv ``` ### 2.6.1 Subfolder `1.diversification.BAMM/` **Task**: Modeling speciation and extinction rate dynamics using BAMM with estimated species richness and currently known species richness. Genus-specific sampling fractions are applied. #### 2.6.1.1 Folder `May-14-2024.Estimated/2.partitioned/treePL/` **Task**: Modeling speciation and extinction rate dynamics using BAMM with estimated species richness. The analyses are run on three selected six-gene trees: ``` . ├── ClipKIT_gappy/ ├── ClipKIT_kpic-smart-gap/ └── ClipKIT_smart-gap/ ``` For instance, in the `ClipKIT_smart-gap/` folder, there are the following folders: ``` . ├── 1.get_prior ├── 2.run1 └── 3.run2 ``` We calculate the priors of related parameters (`lambdaInitPrior`, `lambdaShiftPrior`, and `muInitPrior` ) in the `1.get_prior/` folder: ```bash $ cat get_prior.R library(BAMMtools) setBAMMpriors(read.tree("clean.input.ultrametric.tree")) ``` And then run the analysis in `2.run1/` and `3.run2/` after modify the parameters in the corresponding `control_diversification.txt` file: ```bash $ cd ../2.run1 $ cat b.sh /home/gmeng/soft/bamm-v2.5.0/build/bamm -c control_diversification.txt ``` Note that, in the `control_diversification.txt` file, we use: ``` modeltype = speciationextinction # Specify "speciationextinction" or "trait" analysis treefile = ../1.get_prior/clean.input.ultrametric.tree # File name of the phylogenetic tree to be analyzed sampleProbsFilename = ../1.get_prior/EstimatedNum.tip.group.frac.tsv # File name containing species-specific sampling fractions expectedNumberOfShifts = 50 # prior on the number of shifts in diversification # Suggested values: # expectedNumberOfShifts = 1.0 for small trees (&lt; 500 tips) # expectedNumberOfShifts = 10 or even 50 for large trees (&gt; 5000 tips) ``` The `../1.get_prior/EstimatedNum.tip.group.frac.tsv` file specifies the genus-specific sampling fractions. Now, we have: ``` $ ls *nb step2_Analysis_of_rate_shifts.ipynb step3_rate_through_time_by_groups.ipynb ``` In the Jupyter Notebook file `step2_Analysis_of_rate_shifts.ipynb`, we check if the MCMC chains converge, draw mean phylorate plots (i.e., speciation, extinction, net diversification rates), calculate clade-specific evolutionary rates, etc., following [http://bamm-project.org/documentation.html](http://bamm-project.org/documentation.html). Furthermore, we draw rate-through-time plots in the Jupyter Notebook file `step3_rate_through_time_by_groups.ipynb`. #### 2.6.1.2 Folder `May-14-2024.Current/2.partitioned/treePL/` **Task**: Modeling speciation and extinction rate dynamics using BAMM with currently known species richness. Similar to `May-14-2024.Estimated/2.partitioned/treePL/`. ### 2.6.2 Subfolder `2.trait_evolution.BAMM.subfam-level-sampling-fraction/` **Task**: Similar to `1.diversification.BAMM/`, but with subfamily-specific sampling fractions. ### 2.6.3 Subfolder `3.model_fitting/` **Task**: Perform model fitting analyses for the four traits. ``` . ├── 1.body_size/ ├── 2.leg_length/ ├── 3.leg_length.to.body_size.ratio/ └── 5.met_to_tibia.ratio/ ``` For instance, for the relative leg length (`3.leg_length.to.body_size.ratio/`), we run ```bash Rscript -e "rmarkdown::render('trait_evo_model_fitting.Rmd')" ``` In the `trait_evo_model_fitting.Rmd` file, we test multiple trait evolution models and find the best-fitted model: ```R EvoModels &lt;- c("BM","OU","EB","rate_trend","lambda","kappa","delta","mean_trend","white") fitEvoModels &lt;- lapply(EvoModels, FUN = function(x) fitContinuous(tree, trait_vector, SE=NA, model = x, ncores = 8, control = list(niter=30000, method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent", "subplex")))) AICs &lt;- sapply(fitEvoModels, FUN = function(x) x$opt$aicc) names(AICs) &lt;- EvoModels AICtable &lt;- aicw(AICs) AICtable ``` ### 2.6.4 Subfolder `4.phylogenetic_signal/` **Task**: Calculating phylogenetic signals of the four traits with Lambda and Blomberg’s K. For instance, check the Jupyter Notebook file: `3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/phylogenetic_signal.ipynb` ### 2.6.5 Subfolder `5.ancestral_state/` **Task**: Perform ancestral state reconstruction for the four traits. ``` . ├── 1.body_size ├── 2.leg_length ├── 3.leg_length.to.body_size.ratio └── 5.met_to_tibia.ratio ``` Let's use relative leg length and the ClipKIT_smart-gap tree as an example (`5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/`):: The following R scripts are used for constructing the ancestral states of the trait based on different trait evolution models (BM, EB, OU, trend): ```bash $ ls -l ace.*.r | awk '{print $NF}' ace.BM.r ace.EB.r ace.OU.r ace.trend.r ``` For instance, here we test the BM model: ```R tree_file &lt;- "input.tre" trait_file &lt;- "input.trait.txt" tree &lt;- read.tree(file=tree_file) length(tree$tip.label) df &lt;- read.table(trait_file) head(df) trait_vector0 &lt;- as.vector(df[,2]) head(trait_vector0) names(trait_vector0) &lt;- df[,1] head(trait_vector0) matched_tree_trait_vector &lt;- match.phylo.data(tree, trait_vector0) trait_vector &lt;- matched_tree_trait_vector$data fit &lt;- ace(trait_vector, tree, type="continuous", method="REML", CI=TRUE, model="BM") save(fit, file='fit.BM.Rdata') ``` #### Phenograms After obtaining the ancestral states, we draw the phenograms for all species and subfamilies, respectively. For the subfamily phenogram, we color each subfamily in the whole tree (grey background), so we can see how the trait of each subfamily evolved: ```R # group information specified by the file 'code.subfam.tsv' subfam_group &lt;- read.table('code.subfam.tsv', header = TRUE) subfam_group_2 &lt;- list( Artema=subfam_group[subfam_group$subfam == "Artema", ]$code, Priscula = subfam_group[subfam_group$subfam == "Priscula", ]$code, Arteminae = subfam_group[subfam_group$subfam == "Arteminae", ]$code, Modisiminae = subfam_group[subfam_group$subfam == "Modisiminae", ]$code, Ninetinae = subfam_group[subfam_group$subfam == "Ninetinae", ]$code, Pholcinae = subfam_group[subfam_group$subfam == "Pholcinae", ]$code, Smeringopinae = subfam_group[subfam_group$subfam == "Smeringopinae", ]$code) tree4 &lt;- groupOTU(tree3, subfam_group_2, 'group') # create a 'group' item. pdf(file = "trait.ace.BM.phenogram.color-by-subfam.pdf", width = 8, height = 12) p &lt;- ggtree(tree4, aes(color=group), yscale = "trait", alpha=.9) revts(p) + theme_minimal() dev.off() ``` However, the ancestral state reconstruction step can take a long time to finish, and I want to test different visualization methods. One solution is to first save the ancestral states into a file (e.g., the R code `save(fit, file='fit.BM.Rdata')` in the `ace.BM.r` file). In the end, we use the R script `5.ancestral_state/4.m_tibia_1.to.body_size.ratio/2.partitioned/treePL/ClipKIT_kpic-smart-gap/2nd_phenogram-highlight/2nd_draw_phenogram.highlight.r` to achieve better visualization: **Linear Phenograms** ```bash $ cd 5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/2nd_phenogram-highlight $ Rscript /path/to/5.ancestral_state/4.m_tibia_1.to.body_size.ratio/2.partitioned/treePL/ClipKIT_kpic-smart-gap/2nd_phenogram-highlight/2nd_draw_phenogram.highlight.r -f ../fit.BM.Rdata -p BM ``` **Circular Phenograms** Using the R script `5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/recoloring_phenogram.circular/recoloring_phenogram.circular.R`, we can draw circular phenograms, coloring the tip labels by subfamily : ```bash $ cd 5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/recoloring_phenogram.circular $ Rscript $script_file -f ../fit.BM.Rdata -p BM ``` Furthermore, we also **separate the ancestral states for different subfamilies** using the Jupyter Notebook `5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/trait_std_slop_by_group.ipynb`, which generates a file named `BM.internal-trait_by_group.csv`. ```bash $ head BM.internal-trait_by_group.csv "trait","clade" 9.4,"Pholcinae" 12.7692307692308,"Pholcinae" 14.6923076923077,"Pholcinae" 15.8888888888889,"Pholcinae" 20.2307692307692,"Pholcinae" 19.9166666666667,"Pholcinae" 9.4,"Pholcinae" 7.35849056603773,"Pholcinae" 10.5454545454545,"Pholcinae" ``` #### Trait density plot With the file `BM.internal-trait_by_group.csv`, we can have a look at the overall distribution of the traits. We can do this by drawing a density plot (see `trait_density_plot/draw.trait_density_plot.R`). ```R trait_file &lt;- "../BM.internal-trait_by_group.csv" ylab &lt;- "Clade" xlab &lt;- expression(ratio == frac(leg~length, body~size)) # Customize x-axis title with ratio # Define the colors for each clade (color-blind friendly palette) clade_colors &lt;- c( "Smeringopinae" = "#0072B2", # Blue "Pholcinae" = "#E69F00", # Orange "Arteminae" = "#F0E442", # Yellow "Modisiminae" = "#009E73", # Green "Priscula" = "#D55E00", # Red "Artema" = "#CC79A7", # Pink "Ninetinae" = "#56B4E9", # Sky Blue "Pholcidae" = "#984EA3" # Purple ) ## 3. trait data trait_df &lt;- read.csv(trait_file) head(trait_df) # Set the order of clades trait_df$clade &lt;- factor(trait_df$clade, levels = rev(c("Smeringopinae", "Pholcinae", "Arteminae", "Modisiminae", "Ninetinae", "Artema", "Priscula", "Pholcidae"))) # Specify the desired order p &lt;- ggplot(trait_df, aes(x = trait, y = clade, fill = clade)) + geom_density_ridges(alpha = 0.8) + scale_fill_manual(values = clade_colors) + # Set custom colors theme_ridges() + theme( legend.position = "none", # Remove the legend axis.title.x = element_text(size = 18, hjust = 0.5), # Set x-axis title size and center it axis.title.y = element_text(size = 18, hjust = 0.5), # Set y-axis title size and center it axis.text.x = element_text(size = 17), # Set x-axis text size axis.text.y = element_text(size = 17), # Set y-axis text size ) + labs( x = xlab, # Customize x-axis title y = ylab # Customize y-axis title ) ggsave(filename = "trait_density_plot.pdf", width = 6, height = 5, plot = p) ``` #### Correlations between disparity and diversification rates **Task**: Test for the correlations between trait disparity (measured as the standard variance of the trait values) and speciation/extinction/net diversification rates. The 'rate' files have two versions: (1) Genus-specific sampling fractions are determined by estimated species richness; (2) Genus-specific sampling fractions are determined by currently known species. Please refer to the section **2.6.1 Subfolder** **`1.diversification.BAMM/`**. Therefore, we have the following R script files: ```bash $ ls *r | grep trait1-trait2 Current.trait1-trait2-regression.diversification.r Current.trait1-trait2-regression.extinction.r Current.trait1-trait2-regression.speciation.r Estimated.trait1-trait2-regression.diversification.r Estimated.trait1-trait2-regression.extinction.r Estimated.trait1-trait2-regression.speciation.r ``` **Relationship between trait disparity and speciation rate** For instance, in the `Estimated.trait1-trait2-regression.speciation.r` file: ```R tree_file &lt;- 'input.tre.monophyletic_subfamily.collapsed.tre' rate_file &lt;- '../../../../../1.diversification.BAMM/May-14-2024.Estimated/2.partitioned/treePL/ClipKIT_smart-gap/all_subfamilies_df.csv' trait_file &lt;- "BM.internal-trait_by_group.csv" ylab &lt;- "Standardized Standard Variance" xlab &lt;- "Standardized Speciation Rate" # Define the colors for each clade (color-blind friendly palette) clade_colors &lt;- c( "Smeringopinae" = "#0072B2", # Blue "Pholcinae" = "#E69F00", # Orange "Arteminae" = "#F0E442", # Yellow "Modisiminae" = "#009E73", # Green "Priscula" = "#D55E00", # Red "Artema" = "#CC79A7", # Pink "Ninetinae" = "#56B4E9", # Sky Blue "Pholcidae" = "#984EA3" # Purple ) # read the subfamily-level tree subfam_tree &lt;- read.tree(tree_file) # read the rate file rate_df &lt;- read.csv(rate_file) # it has below columns: #'speciation_rate''extinction_rate''clade''diversification_rate' # Calculate the mean speciation rate for each clade rate_df_2 &lt;- aggregate(speciation_rate ~ clade, data = rate_df, mean) # read the trait file trait_df &lt;- read.csv(trait_file) # Calculate the variance of traits within each clade trait_variance_df &lt;- aggregate(trait ~ clade, data = trait_df, var) names(trait_variance_df)[2] &lt;- "trait_variance" trait_variance_df$standard_variance &lt;- sqrt(trait_variance_df$trait_variance) # now 'trait_variance_df' has three columns: # clade trait_variance standard_variance # Combine the trait variance and speciation rate dataframes by 'clade' combined_df &lt;- merge(trait_variance_df, rate_df_2, by = "clade") # prepare trait data for correlation test trait_1 &lt;- combined_df[['standard_variance']] trait_1 &lt;- setNames(trait_1, combined_df[['clade']]) trait_1 spp &lt;- names(trait_1) trait_2 &lt;- combined_df[['speciation_rate']] trait_2 &lt;- setNames(trait_2, combined_df[['clade']]) trait_2 # Standardize variables trait_1 &lt;- scale(trait_1) trait_2 &lt;- scale(trait_2) # PGLS test # BM model corBM &lt;- corBrownian(phy=subfam_tree,form=~spp) pgls.fit.BM &lt;- gls(trait_1~trait_2,correlation=corBM, method='ML') pgls.fit.BM.summary &lt;- summary(pgls.fit.BM) # Extract coefficient and p-value for trait_2 coefficient_trait_2 &lt;- pgls.fit.BM.summary$tTable[trait_2_row_index, "Value"] p_value_trait_2 &lt;- pgls.fit.BM.summary$tTable[trait_2_row_index, "p-value"] ... ## draw BM pdf("Estimated.trait1_trait2.scatter.reg.PGLS.BM.std.speciation.pdf") ## set margins par(mar = c(6.1, 6.1, 1.1, 1.1), mgp = c(4, 1, 0)) ## graph scatterplot of contrasts plot(trait_1~trait_2, ylab=ylab, xlab=xlab, pch=21, col="black", # Outline color bg=clade_colors[spp], # Fill color based on clades cex=2, las=1, cex.axis=1.5, cex.lab=2, bty="n") text(trait_2, trait_1, labels = spp, cex = 1.5, pos = text_positions, xpd = FALSE) # Add text for coefficient_trait_2 and p_value_trait_2 text(x = 0.02, y = -0.1, cex = 1.5, labels = paste("Slop :", round(coefficient_trait_2, 3), ", p-value:", round(p_value_trait_2, 3)), pos = 4, col = "#984EA3") ## add gridlines to the plot abline(h=0,lty="dotted") abline(v=0,lty="dotted") ## reset graphing limits of the plot to the ## x/y range of our PICs clip(min(trait_2),max(trait_2), min(trait_1),max(trait_1)) # clip(x_min, x_max, y_min, y_max) ## graph our fitted line # a, b : single values specifying the intercept and the slope of the line abline(a = coef(pgls.fit.BM)[1], b = coef(pgls.fit.BM)[2],lwd=2,col="#984EA3") dev.off() ``` A resulting file `Estimated.trait1_trait2.scatter.reg.PGLS.BM.std.speciation.pdf` will be generated, showing the relationship between speciation rate and standard variance ('disparity'). ### 2.6.6 Subfolder `6.relationships/` **Task**: Using the QuaSSE/ES-Sim/PGLS method to test for the correlations between trait evolution and species diversification, among others. This folder contains the following subfolders: ``` . ├── 2.trait-trait.correlation └── 3.trait-dependent-diversification ``` #### 2.6.6.1 `2.trait-trait.correlation/` **Task**: Test for correlations between different traits using PGLS, and get trait distribution over microhabitats (leaf, space, and ground) ``` . ├── 1.body_size.vs.leg_length ├── 2.body_size.vs.leg_length.to.body_size.ratio ├── 3.body_size.vs.m_tibia_1.to.body_size.ratio ├── 4.leg_length.vs.m_tibia_1.to.body_size.ratio ├── 5.body_size.vs.raw-microhabitat ├── 6.body_size.vs.ratio-microhabitat ├── 7.leg_length.vs.raw-microhabitat ├── 8.leg_length.vs.ratio-microhabitat ├── 9.leg_length.to.body_size.ratio.vs.raw-microhabitat ├── 10.leg_length.to.body_size.ratio.vs.ratio-microhabitat ├── 11.m_tibia_1.to.body_size.ratio.vs.raw-microhabitat ├── 12.m_tibia_1.to.body_size.ratio.vs.ratio-microhabitat ├── 13.raw-microhabitat.vs.ratio-microhabitat └── 14.species-richness_vs_trait_std ``` Here are several examples. **Allometric relationship between body size and leg length** ```bash $ cd 6.relationships/2.trait-trait.correlation/1.two-traits/1.body_size.vs.leg_length/2.partitioned/treePL/ClipKIT_smart-gap ``` The script file `two-traits-corelation.R`: ```R treefile &lt;- "input.tre" trait_file &lt;- 'clean.trait1_trait2.txt' ylab &lt;- "log2(body size)" xlab &lt;- "log2(leg length)" tree &lt;- read.tree(file=treefile) length(tree$tip.label) tree trait_df &lt;- read.table(trait_file) head(trait_df) # prepare trait data trait_1 &lt;- trait_df[['V2']] trait_1 &lt;- setNames(trait_1, trait_df[['V1']]) head(trait_1) trait_2 &lt;- trait_df[['V3']] trait_2 &lt;- setNames(trait_2, trait_df[['V1']]) head(trait_2) # ... # ... spp &lt;- names(trait_1) # BM model corBM &lt;- corBrownian(phy=tree,form=~spp) corBM pgls.fit.BM &lt;- gls(trait_1~trait_2,correlation=corBM, method='ML') pgls.fit.BM.summary &lt;- summary(pgls.fit.BM) pgls.fit.BM.summary coef(pgls.fit.BM) abs(coef(pic.fit)[1]-coef(pgls.fit.BM)[2]) summary(pic.fit)$coefficients[1,4] summary(pgls.fit.BM)$tTable[2,4] coef(pgls.fit.BM) ## draw BM pdf("trait1_trait2.scatter.reg.PGLS.BM.pdf") ## set margins par(mar=c(5.1,5.1,1.1,1.1)) ## graph scatterplot of contrasts plot(trait_1~trait_2, ylab=ylab, xlab=xlab, pch=21,bg="gray",cex=1.2,las=1, cex.axis=0.7,cex.lab=0.9,bty="n") ## add gridlines to the plot abline(h=0,lty="dotted") abline(v=0,lty="dotted") ## reset graphing limits of the plot to the ## x/y range of our PICs clip(min(trait_2),max(trait_2), min(trait_1),max(trait_1)) ## graph our fitted line # a, b : single values specifying the intercept and the slope of the line abline(a = coef(pgls.fit.BM)[1], b = coef(pgls.fit.BM)[2],lwd=2,col="red") dev.off() ``` For better visualization, we can color the data points by subfamily, which can be done with the script `6.relationships/2.trait-trait.correlation/1.two-traits/1.body_size.vs.leg_length/2.partitioned/treePL/ClipKIT_smart-gap/coloring_data_points/PGLS_by_subfam.v2.1.R`. **Relationship between trait disparity and species richness** ```bash $ cd 6.relationships/2.trait-trait.correlation/1.two-traits/14.species-richness_vs_trait_std/3.leg_length.to.body_size.ratio/ClipKIT_smart-gap/ ``` We test the relationship with both estimated species richness and currently known species richness: ``` . ├── 1.Estimate └── 2.Current-And-Undescribed ``` In the file, for instance, we have: ```R # smart-gap tree_file &lt;- "/path/to/5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/input.tre.monophyletic_subfamily.collapsed.tre" trait_1_file &lt;- "/path/to/5.ancestral_state/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap/BM.internal-trait_by_group.csv" trait_2_file &lt;- "../Estimated.subfam.count.tsv" ylab &lt;- "Standardized Standard Variance" xlab &lt;- "Standardized Species richness" # Define the colors for each clade (color-blind friendly palette) clade_colors &lt;- c( "Smeringopinae" = "#0072B2", # Blue "Pholcinae" = "#E69F00", # Orange "Arteminae" = "#F0E442", # Yellow "Modisiminae" = "#009E73", # Green "Priscula" = "#D55E00", # Red "Artema" = "#CC79A7", # Pink "Ninetinae" = "#56B4E9", # Sky Blue "Pholcidae" = "#984EA3" # Purple ) ## 1. subfamily-level tree subfam_tree &lt;- read.tree(tree_file) subfam_tree str(subfam_tree) ## 3. trait data trait_1_df &lt;- read.csv(trait_1_file) head(trait_1_df) ## 3.1 Calculate the variance of traits within each clade trait_1_variance_df &lt;- aggregate(trait ~ clade, data = trait_1_df, var) names(trait_1_variance_df)[2] &lt;- "trait_variance" trait_1_variance_df$standard_variance &lt;- sqrt(trait_1_variance_df$trait_variance) trait_1_variance_df # now it has three columns: # clade trait_variance standard_variance trait_2_df &lt;- read.table(trait_2_file, col.names = c("clade", "species_num")) trait_2_df combined_df &lt;- merge(trait_1_variance_df, trait_2_df, by = "clade") combined_df ### prepare trait data trait_1 &lt;- combined_df[['standard_variance']] trait_1 &lt;- setNames(trait_1, combined_df[['clade']]) trait_1 spp &lt;- names(trait_1) trait_2 &lt;- combined_df[['species_num']] trait_2 &lt;- setNames(trait_2, trait_2_df[['clade']]) trait_2 # Standardize variables trait_1 &lt;- scale(trait_1) trait_2 &lt;- scale(trait_2) # BM model corBM &lt;- corBrownian(phy=subfam_tree,form=~spp) corBM pgls.fit.BM &lt;- gls(trait_1~trait_2,correlation=corBM, method='ML') pgls.fit.BM.summary &lt;- summary(pgls.fit.BM) pgls.fit.BM.summary # Find the row index for "trait_2" trait_2_row_index &lt;- which(rownames(pgls.fit.BM.summary$tTable) == "trait_2") # Extract coefficient and p-value for trait_2 coefficient_trait_2 &lt;- pgls.fit.BM.summary$tTable[trait_2_row_index, "Value"] coefficient_trait_2 p_value_trait_2 &lt;- pgls.fit.BM.summary$tTable[trait_2_row_index, "p-value"] p_value_trait_2 # Example vector specifying text positions for each point # 1 = below, 2 = left, 3 = above, 4 = right text_positions &lt;- rep(4, length(trait_1)) # Initialize all positions to 4 (right) # Adjust positions as needed based on data # Here you can customize this vector based on specific conditions for (i in 1:length(trait_1)) { if (spp[i] == "Modisiminae" || spp[i] == "Pholcinae") { text_positions[i] = 2 } } ## draw BM pdf("Estimated.trait1_trait2.scatter.reg.PGLS.BM.std.pdf") ## set margins par(mar = c(6.1, 6.1, 1.1, 1.1), mgp = c(4, 1, 0)) ## graph scatterplot of contrasts plot(trait_1~trait_2, ylab=ylab, xlab=xlab, pch=21, col="black", # Outline color bg=clade_colors[spp], # Fill color based on clades cex=2, las=1, cex.axis=1.5, cex.lab=2, bty="n") text(trait_2, trait_1, labels = spp, cex = 1.5, pos = text_positions, xpd = FALSE) # Add text for coefficient_trait_2 and p_value_trait_2 text(x = 0.02, y = -0.1, cex = 1.5, labels = paste("Slop :", round(coefficient_trait_2, 3), ", p-value:", round(p_value_trait_2, 3)), pos = 4, col = "#984EA3") ## add gridlines to the plot abline(h=0,lty="dotted") abline(v=0,lty="dotted") ## reset graphing limits of the plot to the ## x/y range of our PICs clip(min(trait_2),max(trait_2), min(trait_1),max(trait_1)) # clip(x_min, x_max, y_min, y_max) ## graph our fitted line # a, b : single values specifying the intercept and the slope of the line abline(a = coef(pgls.fit.BM)[1], b = coef(pgls.fit.BM)[2],lwd=2,col="#984EA3") dev.off() ``` **Leg length distribution over microhabitats**: ```bash $ cd 6.relationships/2.trait-trait.correlation/1.two-traits/7.leg_length.vs.raw-microhabitat/2.partitioned/treePL/ClipKIT_smart-gap ``` ```R pdf("trait-distribution-by-microhabitat.pdf") cols &lt;- c("#F76D5E", "#FFFFBF", "#72D8FF") ggplot() + geom_density(data=trait_df, aes(x = trait_1), colour = "red", linetype = "dashed") + geom_density(data=trait_df, aes(x = trait_1, fill = microhabitat), alpha = 0.6) + scale_fill_manual(values = cols) + xlab(trait_1_name) dev.off() ``` #### 2.6.6.2 `3.trait-dependent-diversification/` **Task**: Test for correlations between trait evolution and species diversification. I will explain this with relative leg length-dependent diversification as an example: ```bash $ cd 6.relationships/3.trait-dependent-diversification/1.one-trait-dependent/5.leg_length.to.body_size.ratio.vs.divrsification/2.partitioned/treePL/ClipKIT_smart-gap ``` **QuaSSE** The analysis here basically follows the following tutorial * *Analysing diversification with diversitree* (Rich FitzJohn, with GeoSSE by Emma Goldberg; [https://www.zoology.ubc.ca/prog/diversitree/doc/diversitree-tutorial.pdf](https://www.zoology.ubc.ca/prog/diversitree/doc/diversitree-tutorial.pdf) The QuaSSE (Quantitative State Speciation and Extinction) model in the R package diversitree (v0.9-20) (FitzJohn 2010), which simultaneously analyzes the relationship between changes in continuous traits and diversification and can suffer from model misspecification and phylogenetic pseudoreplication (Harvey and Rabosky 2018). A global sampling fraction was applied, since the model does not support clade-specific sampling fractions. We fitted the speciation and extinction rates as constant (`1.nodrift.mle.c.R`), linear (`*.mle.l.R`), sigmoidal (`*.mle.s.R`), and hump-shaped (`*.mle.h.R`) functions of the traits, respectively. We also tested the models without the drift parameter (i.e., assuming a constant mean of trait values through time). Finally, we performed a likelihood ratio test to compare the model fits. ```bash $ ls -l *R | awk '{print $NF}' 1.nodrift.mle.c.R 2.1.mle.s.R 2.2.mle.l.R 2.3.mle.h.R 3.1.mle.s.R 3.2.mle.l.R 3.3.mle.h.R 4.1.mle.s.R 4.2.mle.l.R 4.3.mle.h.R ``` The `1.nodrift.mle.c.R` assumes speciation ('lambda') and extinction ('mu') rates are a constant function of trait (`f.c &lt;- make.spiders(constant.x, constant.x)`). The `2.1.mle.s.R` assumes the speciation rate is a sigmoidal function of a trait, while the extinction rate is a constant function of a trait (`f.s &lt;- make.spiders(sigmoid.x, constant.x)`), combined with the presence or absence of the drift parameter: ```R mle.s &lt;- find.mle(nodrift(f.s), p.s, control=control, verbose=0) mle.d.s &lt;- find.mle(f.s, coef(mle.s, TRUE), control=control, verbose=0) ``` The `3.1.mle.s.R` assumes the extinction rate is a sigmoidal function of a trait, while the speciation rate is a constant function of a trait (`f2.s &lt;- make.spiders(constant.x, sigmoid.x)`), combined with the presence or absence of the drift parameter: ```R mle2.s &lt;- find.mle(nodrift(f2.s), p2.s, control=control, verbose=0) mle2.d.s &lt;- find.mle(f2.s, coef(mle2.s, TRUE), control=control, verbose=0) ``` The `4.1.mle.s.R` assumes both speciation rate and extinction rate are a sigmoidal function of a trait (`f3.s &lt;- make.spiders(sigmoid.x, sigmoid.x)`), combined with the presence or absence of the drift parameter: ```R mle3.s &lt;- find.mle(nodrift(f3.s), p3.s, control=control, verbose=0) mle3.d.s &lt;- find.mle(f3.s, coef(mle3.s, TRUE), control=control, verbose=0) ``` Note that the estimated parameter values from a previous step serve as the initial values of the subsequent steps (i.e., `1.nodrift.mle.c.R` -&gt; `2.1.mle.s.R` -&gt; `3.1.mle.s.R` -&gt; `4.1.mle.s.R`). Each step can take a very long time to run. This is why I save the results from the previous step into an Rdata file (e.g., `save.image(file='until.nodrift.mle.c.Rdata')`) and load them again at the subsequent step (e.g., `load('until.nodrift.mle.c.Rdata')`). Finally, model comparison is performed with the `anova()` function in R. Please refer to the Jupyter Notebook `single_partition.annova.ipynb` for details. **ES-Sim test** The tip-rate correlation ES-sim test, which does not require a fully parameterized model of diversification and trait evolution, and directly examines the correlation between trait variation across the tips of a phylogeny and tip-specific speciation rate estimates, and is thus believed to have rare false inferences (Harvey and Rabosky 2018). We used 20,000 simulations to build the null distribution of trait-speciation associations for significance testing: ```R es_result &lt;- essim(tree, trait_vector, nsim = 20000) ``` Please refer to [https://github.com/mgharvey/ES-sim](https://github.com/mgharvey/ES-sim) for more details. ### 2.6.7 Subfolder `7.disparity.28-05-2024/` **Task**: Draw disparity through time plots with the function `dtt()` for each subfamily and each trait. For instance, ```bash $ cd 7.disparity.28-05-2024/3.leg_length.to.body_size.ratio/2.partitioned/treePL/ClipKIT_smart-gap Rscript ../../../../subfamily_disparity.r ``` For better visualization, I also save the related data into a file (see the file `subfamily_disparity.r`): ```R save.image(file='my_subfamily_disparity.RData') ``` With the data file `my_subfamily_disparity.RData`, we redraw the disparity through time plot using the Jupyter Notebook file `redraw-subfam-disparity.ipynb`. The resulting files are `all_dtt_plots.colored.pdf` and `all_dtt_plots.r2.colored.pdf`, where different subfamilies are colored with different colors. ### 2.6.8 Subfolder `8.linage_through_time/` Task: Draw a Lineage through time plot for the Pholcidae family using the function `ltt.plot()`. For instance, use the Jupyter Notebook file `6_comparative_analysis/8.linage_through_time/2.partitioned/treePL/ClipKIT_smart-gap/llt_plot.ipynb`. ### 2.7 Folder `7_UCE_and_six-gene_phylogeny/` **Task**: Perform phylogenetic analysis with UCE+Sanger data. **Result file**: `combined_UCE+Sanger_trees.pdf` **Main conclusion**: The phylogenetic trees from UCE+Sanger data are highly consistent with our six-gene trees, and all misplacements of taxa are found to be represented by CO1 only; thus, our conclusions on the correlations among evolutionary flexibility and other variables are robust. In this section, we performed: 1. **unpartitioned** analysis based on **ALL UCE loci + Sanger data**, using IQ-TREE with the optimal model determined by ModelFinder (Kalyaanamoorthy et al. 2017) in IQ-TREE automatically. 2. **partitioned** anlaysis based on **55 UCE loci + Sanger data**. The 55 UCEs were selected from matrix 3, and they have **the most parsimonious sites**. The tree was constructed with RAxML-NG and the model GTR+R4+I+FO. You can find the commands and alignments, resulting tree files (`*.treefile`) at the corresponding folders. ``` ./ ├── 1_unpartitioned │   ├── matrix_3.gappy │   │   ├── FcC_supermatrix.phy │   │   ├── iqtree.sh │   │   ├── m3_unpartitioned_gappy.longName.tre │   │   └── m3_unpartitioned_gappy.treefile │   ├── matrix_3.kpic-smart-gap │   │   ├── FcC_supermatrix.phy │   │   ├── iqtree.sh │   │   ├── m3_unpartitioned_kpic-smart-gap.longName.tre │   │   └── m3_unpartitioned_kpic-smart-gap.treefile │   └── matrix_3.smart-gap │   ├── FcC_supermatrix.phy │   ├── iqtree.sh │   ├── m3_unpartitioned_smart-gap.longName.tre │   └── m3_unpartitioned_smart-gap.treefile ├── 2_partitioned │   ├── matrix_3.55UCEs.gappy │   │   ├── 55UCE-Sanger.partitioned.RAxML-NG.gappy.longName.tre │   │   ├── 55UCE-Sanger.partitioned.RAxML-NG.gappy.tree │   │   ├── merged.m3.six-gene.FcC_supermatrix_partition.txt │   │   ├── merged.m3.six-gene.FcC_supermatrix.phy │   │   └── raxml.sh │   ├── matrix_3.55UCEs.kpic-smart-gap │   │   ├── 55UCE-Sanger.partitioned.RAxML-NG.kpic-smart-gap.longName.tre │   │   ├── 55UCE-Sanger.partitioned.RAxML-NG.kpic-smart-gap.tree │   │   ├── merged.m3.six-gene.FcC_supermatrix_partition.txt │   │   ├── merged.m3.six-gene.FcC_supermatrix.phy │   │   └── raxml.sh │   └── matrix_3.55UCEs.smart-gap │   ├── 55UCE-Sanger.partitioned.RAxML-NG.smart-gap.longName.tre │   ├── 55UCE-Sanger.partitioned.RAxML-NG.smart-gap.tree │   ├── merged.m3.six-gene.FcC_supermatrix_partition.txt │   ├── merged.m3.six-gene.FcC_supermatrix.phy │   └── raxml.sh ├── April-9-2024.code.subfam.tsv ├── combined_UCE+Sanger_trees.pdf └── combined.code.subfam.tsv ``` * `April-9-2024.code.subfam.tsv`: relationships between the sample codes and long names * `combined.code.subfam.tsv`: indicating to which subfamily a sample code belongs.