Continuous trait evolution

Data

We switch data for this section, to use a larger network with more taxa: on Xiphophorus fishes from Bastide et al. (2018), available on dryad. More specifically, download

• the time calibrated network
• the trait data on sword index and female preference for a sword, originally from Cui et al. (2013).

The files are also in data_results/cont-traits.

Preparation

If not done already, load the packages needed for this analysis:

using PhyloNetworks
using PhyloPlots
using RCall
using CSV
using DataFrames
using StatsModels

then read in the networks and the trait data:

topologies = readMultiTopology("xiphophorus_networks_calibrated.tre");
net3 = topologies[3]; # we will use the network with 3 reticulations
plot(net3, :R)        # just to look at the network: topology only
plot(net3, :R, useEdgeLength=true) # topology + branch lengths
taxa_net = tipLabels(net3) # extract list of taxa

Next, we will delete the rows in the trait data, for the taxa that are missing in the tree.

dat = CSV.read("xiphophorus_morphology_Cui_etal_2013.csv", copycols=true)
for i in reverse(1:nrow(dat))
    j = findfirst(isequal(dat[:tipNames][i]), taxa_net)
    if j==nothing # taxon not found in network
        println("taxon ", dat[:tipNames][i], " (row $i) not found in network")
        deleterows!(dat, i) # error
    end
end

The snippet above should work fairly generically. Here, it tells us that "taxon Xnezahualcoyotl (row 19) not found in network".

Ancestral state prediction

sword index:

dat_si = dat[:, [:sword_index, :tipNames]]
AS_si = ancestralStateReconstruction(dat_si, net3)
R"par(mar=c(.5,.5,.5,.5))";
plot(net3, :R, nodeLabel=expectationsPlot(AS_si), xlim=[0,25]);
plot(net3, :R, nodeLabel=predintPlot(AS_si), useEdgeLength=true, xlim=[-3,26]);
# compare ancestral state prediction interval to:
extrema(dat[:sword_index])

female preference for a sword:

dat_fp = dat[:, [:preference, :tipNames]]
AS_fp = ancestralStateReconstruction(dat_fp, net3)
plot(net3, :R, nodeLabel=expectationsPlot(AS_fp), xlim=[0,25]);
plot(net3, :R, nodeLabel=predintPlot(AS_fp), useEdgeLength=true, xlim=[-3,26]);
extrema(skipmissing(dat[:preference]))

Phylogenetic signal: Pagel's lambda

lambda_si = phyloNetworklm(@formula(sword_index ~ 1), dat, net3, model="lambda")
lambda_fp = phyloNetworklm(@formula(preference ~ 1),  dat, net3, model="lambda")

Phylogenetic regression / ANOVA

Regression here: does preference influence sword index?

@rlibrary ggplot2 # assumes R package ggplot2 already downloaded & installed
ggplot(dat, aes(x=:preference, y=:sword_index)) + geom_point(alpha=0.9) +
  xlab("female preference") + ylab("sword index");
fit_BM = phyloNetworklm(@formula(sword_index ~ preference), dat, net3)
fit_λ  = phyloNetworklm(@formula(sword_index ~ preference), dat, net3, model="lambda")

Transgressive evolution

regNet = regressorHybrid(net3)
plot(net3, :R, showEdgeNumber=true, xlim=[0,18]);
dat3 = join(dat, regNet, on = :tipNames)
fit0 = phyloNetworklm(@formula(sword_index ~ 1),   dat3, net3)
fit1 = phyloNetworklm(@formula(sword_index ~ sum), dat3, net3)
fit2 = phyloNetworklm(@formula(sword_index ~ shift_24 + shift_37 + shift_45), dat3, net3)
ftest(fit2, fit1, fit0)