Morphological Diversification and Rates of Evolution

This tutorial will show you how to run analyses in Brownie, a package designed by Brian O’Meara (http://www.brianomeara.info/) to estimate and compare rates and means of character evolution in different parts of a phylogeny (see O’Meara et al. 2006).

Installing Brownie

---

Go to here and download BrownieInstaller2.1p2.dmg. This exercise assumes you are running Brownie on Mac OSX with an Intel processor. Using Brownie on a PC is possible, but can be difficult. You will need to install cygwin, install the gsl, and then compile Brownie from source. See the Brownie manual for more details.

Tutorial: Rates of Morphological Evolution

---

Rates of character evolution are not only interesting in their own right but are also a useful tool for other comparative analyses such as ancestral state reconstruction and, in a morphological context, Brownian rates of evolution can be used as phylogenetically corrected estimates of morphological variance (see O'Meara et al. 2004 and Hutcheon & Garland 2004). For more information concerning the comparison of means see discussion of the Ornstein-Uhlenbeck model by Hansen 1997 and Butler and King 2004 as well as the package OUCH in ‘R’. It is also worth noting that code written for ‘R’ by Gavin Thomas and colleagues allows the fitting of models under Brownian motion that allow both the rates and means to vary (see Thomas et al. 2009).

The original Brownie function was a censored rate test under Brownian motion, which allows the comparison of the rates of evolution across two monophyletic sister-clades or between a monophyletic and paraphyletic clade. The censored rate test allows you answer questions such as “Are rates of floral evolution in the large parasitic Rafflesia faster than their smaller relatives?” (Barkman et al. 2008). The newer functions implement the non-censored rate test, which allows you to compare the rates (Brownian motion model) or means (Ornstein-Uhlenbeck model) of evolution in different character states that are mapped (or painted) onto the phylogeny. The non-censored rate test allows you to answer questions such as “Does piscivory limit morphological diversification?” (Collar et al. 2009).

We are using a dataset on parrotfish jaw morphology parrot.nex to answer two questions:

1) Using the non-censored rate test "Does coral excavating/scraping limit morphological evolution in parrotfishes?


2) Using the censored-rate test "Is the evolution of an intramandibular joint in a sub-clade of parrotfishes linked to a increase in the rate of jaw evolution?"

See the papers cited in this tutorial, Brian O’Meara’s website and the Brownie manual http://www.brianomeara.info/brownie/manual for more details.

PLEASE NOTE it is always a good idea to run the analyzes several times even if you don't get warning messages just to make sure that the results are robust as you can sometimes be surprised!

1. Prior to using Brownie

---

The tree and data being input into Brownie need to be in a single nexus file with separate taxa, characters, and trees blocks, as well an optional assumptions block. Since it's unlikely that you will be starting with this type of nexus file, the easiest way to get one ready for analysis with Brownie is to use Mesquite. In Mesquite, you can import data from an EXCEL file containing your character data, link this to a tree file (numerous types of tree formats can be imported into Mesquite), and save this as a nexus that contains most of the information required by Brownie. Additional information on how to generate such a nexus file can be found [Basic ancestral reconstruction with parsimony & maximum likelihood in Mesquite here]. If you have multiple trees, say from a BEAST run, I would recommend you only use 1 tree to generate your Mesquite nexus file as it can take a long time to open and may crash, you can cut and paste the rest of the trees in to the nexus file later. Although there may be some way to add a taxset to an assumptions block in Mesquite this can also be done by manually manipulating the assumptions block in the nexus file output by Mesquite. Today we don’t have to worry about any of this as we are starting with a nexus formatted dataset parrot.nex.

2. Getting started with Brownie

---

a. Open Brownie.

i) If you downloaded the .dmg file, double click on it and install Brownie, then double click on the Brownie icon and choose your input file (you can ignore section 1. b and move straight to 1.c).

ii) If you have the command line version open Terminal (Mac) and navigate to the Brownie folder and move the nexus file with tree and data into that folder. The simplest way to do this may be to simply drag and drop the nexus files into the Brownie before opening the Terminal. Movement of files can also be done in the Terminal. For example, if you want to move the file parrot.nex from your desktop into a folder called Brownie that is also on your Desktop you would type:

>cd Desktop 
>mv parrot.nex Brownie/
>cd Brownie

b. Once you've navigated to the folder containing Brownie you're ready to start the program and execute the example file called parrot.nex. This is done using the command brownie followed by the name of the file you'd like to run.
>brownie parrot.nex

If this didn't work you may also need to add ./ in front of the Brownie instruction (./brownie parrot.nex).
Alternatively you downloaded the installer .dmg you can open Brownie by typing “brownie” and then executing the file once Brownie has started up by typing “exe parrot.nex” or “execute parrot.nex”.

c. Check tree and data (always good practice to see if the data and tree has been read-in correctly etc.)

>showtree 
>blocks

d. Start a log file to store all of the output printed to your screen.

>log file=parrotlog.log

e. Start echoing output to store all commands in a batch file. If you do this, you'll have a stored record of your commands that can subsequently be used to write a Brownie block in the nexus file. In this manner it will be possible to run your analyses automatically when you execute your nexus file.

>echo file=parrotcommand.txt

Below is what you should see on the screen when Brownie starts up and you run the commands on the previous page – if you don’t there is a problem!

3. Non-censored rate test

---

Using stochastic character mapping (see Nielsen 2002, Huelsenbeck et al. 2003, which is amongst other things, a Bayesian approach to ancestral state reconstruction that can take into account uncertainty in phylogenetic reconstruction, branch length estimation and in character reconstruction) or other ancestral state reconstruction methods you can compare the rates (Brownian motion model) or means (Ornstein-Uhlenbeck model) of evolution in different character states, for example “does jaw morphology evolve faster in freshwater or marine fishes?”

Brownie uses the SIMMAP output format (http://www.simmap.com, which implements stochastic character mapping) as its input format to load the ancestral state reconstruction(s), thus it is easiest to do your reconstructions in SIMMAP and then replace the tree(s) in your nexus file with your set of SIMMAP trees. Note that you can also use Brownie’s built-in discrete character optimization to create a set of reconstructions (see http://www.brianomeara.info/brownie/manual#Discretecharacterevolution) but it is relatively untested.

Today we don’t have to worry about any of this as we are starting with a nexus formatted dataset that includes 10 SIMMAP trees: parrot.nex. Please note this is a very small number of SIMMAP trees, as we want the analyses to run quickly, usually you will want 500+ to ensure correct representation of the likely character mappings.

a. Specifying the model

Before executing the non-censored rate test you first need to specify the model you want to use using the ‘model’ command, the default is a single-rate Brownian model (BM1). To see the model options type:

>model ?

BM1 = Brownian motion, one rate parameter
BMS = Brownian motion, with different rate parameters for each state on a tree
BMC = Brownian motion, with one rate parameter for branches with state changes and another for branches without changes
OUSM = Ornstein-Uhlenbeck with one mean per discrete state (attraction and rate parameters constant across tree)
OUCM = Ornstein-Uhlenbeck with independent means on branches with and without changes in a discrete character (attraction and rate parameters constant across tree)

Today we will be using BM1 and BMS, start with BM1 (the default) to answer the question "Does coral excavating limit morphological diversity in parrotfishes?"

>model type=BM1

b. Running a non-censored rate test

To test for different rates of morphological evolution in different parts of the tree, we need to become familiar with Brownie's cont command. To see the options associated with this command type cont ? while Brownie is open (be careful to put a space between the command and the ?)

>cont ?

Taxset: This command is specified in the optional assumptions block of the nexus file. One way to generate this set is to manually edit a nexus file created by Mesquite by specifying an assumption block and a taxset. To see an example of the results of such editing, scroll to the bottom of the parrot.nex file in a text editor. It specifies which taxa in the tree should be used in the analyses, usually when running a non-censored test it is all the species in the tree but it can be used to remove species with missing data etc. For the non-censored rate analysis use the taxset=all, you will also notice an additional taxset labeled intrajoint for use when running the censored-rate test.

Treeloop: This command is used to two Brownie if it should loop across trees (with or without weights) in a nexus file. This looping is used to account for uncertainty in tree topology and branch lengths, or even to weight the average when the output is summarized across input trees with bootstrap support or posterior probabilities.. If you want to use a specific tree use the ‘choose tree’ command to pick the one you want, for example to use the 50th tree in your nexus file type choose tree=50

Charloop: This command can be used to loop across the characters specified in the character block.

File: Specifies a file name in which the output should be saved.

Append and Replace refer to whether the output file should replace or be appended to an already existing file.

**Now we are ready to run our analyses**
We will compare a model that fits a single-rate model to one that allows the rates to vary across the mapped character you need to run two tests: firstly model type=BM1 (1-rate model) which we have already chosen and is the default option and model type=BMS (multiple-rate model). In this example we are comparing the rates of morphological evolution across feeding modes in parrotfishes (0 browser, 1 scraper/excavator), to answer the question Does coral excavating/scraping limit morphological evolution in parrotfishes?

>cont taxset=all treeloop=yes charloop=y file= parrotfish_feedingmodetest.txt 

>model type=BMS

>cont taxset=all treeloop=yes charloop=y file= parrotfish_feedingmodetest.txt append=yes replace=no

>echo stop

>log stop

This should take approx. 90 seconds to run

c. Interpreting the output of a non-censored rate test

Brownie outputs to screen unless you specify an output file in the ratetest command, if you started a log file it will save the screen output. Open the output file parrotfish_feedingmodetest.txt in EXCEL, it will be a bit messy but it is quick to clean it up. Start by naming the column which is named #NAME? to minusLnl, then sort by Model and then character and delete all the rows that are column names at the bottom of the dataset. You should now have 40 lines of data one for each combination of tree (10 trees), model (2 models) and character (2 characters).

We want to compare the single-rate model (BM1) to the equivalent 2-rate model (BMS) and see if one fits the data substantially better, we can do this using a model fitting approach with the Akaike Information Criterion (AIC – see Burnham and Anderson 2002 for an explanation). We calculate the change in the AIC value (∆AIC) for each tree-trait combination (BM1 – BMS) or perhaps better the ∆AICc, which is corrected for small sample sizes. You can then calculate the average rate and ∆AICc across the different feeding mode reconstructions as well as looking at the effect of individual reconstructions (e.g. is there a single reconstruction that leads to the opposite conclusion etc.?). The final EXCEL spreadsheet should look something like this:

From this output we can see that there is no preference for a 1 or 2-rate model in Char 1, which is pc1 (i.e. the first column of data in parrot.nex), as the ∆AICc is <4 (which is the threshold recommended for AIC see Burnham and Anderson 2002) but the single-rate model always gets the lowest AICc score. If you look at the rate estimates you can see a general trend for slower rates in state 0, which is the browsing feeding mode. In Char 2, which is pc2, there is some support for a 2-rate model as the average ∆AICc 3.8 and on average the rate of jaw evolution is x6 slower in browsers (state 0) than excavators/scrapers (state 1).

Does coral excavating/scraping limit morphological evolution in parrotfishes? No, in fact there is an indication that the rate of jaw evolution is faster (and therefore morphological disparity is higher) within excavating/scraping species in pc2.

4. Censored rate test

---

a. Running a censored rate test
This compares the rates of evolution in a monophyletic clade to another monophyletic or paraphyletic clade i.e. rate change only occurs at nodes (in contrast to the cont command which allows changes to occur anywhere along the branch or at the node). It is useful for asking questions about whether a key character has lead to increased rates of evolution. It is ‘censored’ as it removes the branches that join the two clades together as it is not clear which clade to assign these branches to.

To test for different rates of morphological evolution in different parts of the tree, we need to become familiar with Brownie's ratetest command. To see the options associated with this command type ratetest ? while Brownie is open (be careful to put a space between the command and the ?):

Taxset: This command is specified in the optional assumptions block of the nexus file. One way to generate this set is to manually edit a nexus file created by Mesquite by specifying an assumption block and a taxset. To see an example of the results of such editing, scroll to the bottom of the parrot.nex file in a text editor. In the context of the taxset command names one of the two groups of interest and specifies the species names belonging to that group. ratetest will automatically generate another taxset that includes all the species in the tree that are not in the specified taxset. Our example will compare rates among a monophyletic clade to a paraphyletic clade but it can also be used when both clades are monophyletic. In this tutorial we will be using the taxset='intrajoint' (not taxset=all which is used in section 3. Non-censored rate test).

Reps: This command specifies how many times data should be simulated under the null model when doing a parametric bootstrap to assess the significance (p-value) of the likelihood ratio test. This is a good idea if you have small sample sizes but can take a long time for larger datasets. If you are using the model testing approach Reps is unnecessary because your inference will be based on differences in AIC or AICc scores.

Treeloop: This command is used to two Brownie if it should loop across trees (with or without weights) in a nexus file. This looping is used to account for uncertainty in tree topology and branch lengths, or even to weight the average when the output is summarized across input trees with bootstrap support or posterior probabilities.

Charloop: This command can be used to loop across the characters specified in the character block.

File: Specifies a file name in which the output should be saved.

**Now we are ready to run our analyses**

We will compare a model that fits a single-rate to one that allows the rates to vary between the two clades we have specified in taxset. In this example we are comparing the rates of morphological evolution between a clade of parrotfishes that possess an intramandibular joint in their jaws (intrajoint= Scarus, Chlororus, Hipposcarus) to those that don’t (a paraphyletic group), to answer the Is the evolution of an intramandibular joint in a sub-clade of parrotfishes linked to a increase in the rate of jaw evolution?

>ratetest taxset=intrajoint reps=1000 charloop=yes  treeloop=no file=intrajointratetest.txt

c. Interpreting the output of a censored rate test

Brownie outputs to screen unless you specify an output file in the ratetest command, if you started a log file it will save the screen output. Below is the screen output for both characters.

I) This first section of output gives the parameter estimates for each clade in the 2-rate model. Anc state gives the maximum likelihood estimate of the ancestral state at the root node of each clade (Notbaleen or baleen). Rate gives the estimated Brownian motion rate for each clade and –lnL gives the likelihood of that rate.

II) The second section gives the parameter estimates for the 1-rate model. First is –lnL, the likelihood of the single-rate model, followed by two values of interest if you are using the model-fitting approach: the Akaike Information Criterion (AIC) which estimates the goodness of fit of the 1-rate model and AICc which corrects for small sample sizes. It ends with the estimated Brownian motion rate for the single-rate model.

III) This third section summarizes the overall parameters for the 2-rate model. First is –lnL, the likelihood of the two-rate model, which is the sum of the likelihoods of the baleen and Notbaleen Brownian motion models given in the first section. This is followed by the AIC and AICc values for the 2-rate model.

IV) The fourth section summarizes the results and conclusions from the various significance tests and model testing approaches: AIC differences, AICc differences, Χ2 p-value which gives the significance of a likelihood ratio test that compares the likelihood of the 2-rate model with the single-rate model and finally because we ran Reps=1000 the p-value for the parametric bootstrapping of the significance of the likelihood ratio test.

v) Once it has finished looping over the characters a Summary of results is printed to the screen. Brownie uses the notation A,a,B,b for strength of support for a single rate (A’s) or two rate (B’s) model with strong support indicated by CAPITALS.

From this output we can see that there is weak support (parametric bootstrap p-value 0.05 and ∆AICc 1.7) for a 2-rate model in pc1 (char 1): the intramandibular joint clade exhibits a x2.8 faster rate of jaw evolution. The support for a 2-rate model is stronger in pc2 (∆AICc 6.5 and parametric bootstrap p-value 0.005) with the intramandibular joint clade exhibiting an almost 5x faster rate of jaw evolution.

Is the evolution of an intramandibular joint in a sub-clade of parrotfishes linked to a increase in the rate of jaw evolution? Yes, it appears that pc2 and probably pc1 evolve at faster rates within the three genera that possess an intramandibular joint. See presentation by Peter Wainwright for further details on this pattern.

5. Making a batch file

---

This will create a file that will run the above analysis in section 3. when the nexus file is executed – this is particularly useful if submitting it run on a cluster where you are less likely to be able to work interactively with Brownie.

Open in a text editor the commands that we echoed to parrotcommand.txt and open the nexus data file parrot.nex. Copy all the text below #nexus in the parrotcommand.txt file and paste it at the bottom of the text in the parrot.nex file, save this new file as batchparrot.nex and then execute the batchparrot.nex file in Brownie.

>Execute the batchparrot.nex

6. Warnings and errors

---

a. When running the non-censored rate test you may receive the warning:

---

WARNING: Out of 15 optimization starts, 15 were stopped by hitting

the maximum # of iterations. This means that those replicates


may not even have hit the local maximum.

You can increase the maximum number of iterations or decrease the


precision with the NumOpt command. You could also consider


increasing the number of random starts using that same command.

If this happened on a small proportion of replicates, though,


or if the precision (below) is good enough, don't worry about it.

---

Obviously this is not good! However, to see if this is actually affecting your results run the analyses several times and see if the results change substantially (this is what it means by precision), if not you are probably okay. Also, as it states in the warning, you can try increasing the number of random starts and/or iterations using the numopt command e.g.

>numopt randstart=50 iter=5000

>numopt ?

Iter sets the maximum number of iterations of the Nelder-Mead simplex algorithm.

Toler sets the precision of the stopping criterion: what amount
of change in the likelihood is considered small enough to count as zero change.

RandStart sets the number of random starts to use.

StepSize sets the NM step size.

Detail specifies whether or not to have detailed output from numerical optimization

Redo specifies whether to redo reps, which stop due to iteration limits

GiveUpFactor, when redo=yes, is used to tell the software when to stop restarting: when the ratio of unsuccessful to successful starts is > giveupfactor

b. BUS errors, segmentation faults – these are memory errors. The current version of Brownie has some memory leaks so these might be causing Brownie to quit once it has used all the available memory. Try moving to a machine with more memory or increasing your virtual memory!

7. References

---

Barkman, Bendiksby, Lim, Salleh, Nais, Madulid and Schumacher (2008) Accelerated rates of floral evolution at the upper size limit for flowers. Current Biology 18(19) pg 1508-1513.

Burnham, K. P., and D. R. Anderson. 2002. Model Selection and Multimodel
Inference: A practical information theoretic approach. Springer-Verlag, New York.

Butler, M. A., and A. A. King. 2004. Phylogenetic comparative analysis: a
modeling approach for adaptive evolution. American Naturalist 164:683-695.

Collar D. et al. 2009. Piscivory limits diversification of feeding morphology in centrarchid fishes. Evolution 63(3), 1557-1573.

Hansen, T. F. 1997. Stabilizing selection and the comparative analysis of
adaptation. Evolution 51:1341-1351.

Huelsenbeck, J. P., Nielsen, R., Bollback, J. P. 2003. Stochastic mapping of morphological characters. Systematic Biology 52(3), 131-158.

Hutcheon & Garland 2004. Are megabats big? Journal of Mammalian Evolution 11 (3-4), pg 257-277

Nielsen, R. 2002. Mapping mutations on phylogenies. Systematic Biology 51(5), 729-739.

O'Meara B. et al. 2006. Testing for different rates of continuous trait evolution using likelihood. Evolution 60(5), pg 922-933.

Revell, L. 2009. Size-correction and principal components for interspecific comparative studies. Evolution 63(12), 3258-3268.

Thomas, GH, Meiri, S & Phillimore, AB. 2009. Body size diversification in Anolis: novel environments and island effects. Evolution 63(8), 2017-2030.

Wainwright, P. C., D. R. Bellwood, M. W. Westneat, J. R. Grubich, and A. S. Hoey. 2004. A functional morphospace for the skull of labrid fishes: patterns of diversity in a complex biomechanical system. Biological Journal of the Linnean Society 82:1-25.

Comments:

Note: You must be logged in to add comments