diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
index f49313b..0cd3114 100644
--- a/.github/workflows/TagBot.yml
+++ b/.github/workflows/TagBot.yml
@@ -4,6 +4,22 @@ on:
     types:
       - created
   workflow_dispatch:
+    inputs:
+      lookback:
+        default: "3"
+permissions:
+  actions: read
+  checks: read
+  contents: write
+  deployments: read
+  issues: read
+  discussions: read
+  packages: read
+  pages: read
+  pull-requests: read
+  repository-projects: read
+  security-events: read
+  statuses: read
 jobs:
   TagBot:
     if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
diff --git a/Project.toml b/Project.toml
index fff4eb9..bcba59c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PhyloPlots"
 uuid = "c0d5b6db-e3fc-52bc-a87d-1d050989ed3b"
 license = "MIT"
-version = "1.0.1"
+version = "2.0.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -12,9 +12,9 @@ RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
 
 [compat]
 DataFrames = "0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 1.0"
-PhyloNetworks = "0.16, 0.17.0, 1.0.0" # extras to let PN build its documentation of new releases
+PhyloNetworks = "0.17"
 RCall = "0.11, 0.12, 0.13, 0.14"
-julia = "0.7, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7"
+julia = "1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/README.md b/README.md
index cb198c7..1d43361 100644
--- a/README.md
+++ b/README.md
@@ -17,9 +17,11 @@ via [RCall](https://github.com/JuliaInterop/RCall.jl).
 
 ## examples
 
-in PhyloPlots' documentation [manual](https://juliaphylo.github.io/PhyloPlots.jl/stable/man/getting_started/), and  
-in [PhyloNetworks](http://juliaphylo.github.io/PhyloNetworks.jl/latest/)'s documentation, see
-- basic examples
-  [here](http://juliaphylo.github.io/PhyloNetworks.jl/latest/man/snaq_plot/#Network-Visualization-1)
-- examples to annotate edges
-  [here](http://juliaphylo.github.io/PhyloNetworks.jl/latest/man/bootstrap/#support-for-tree-edges)
+- in PhyloPlots' documentation [manual](https://juliaphylo.github.io/PhyloPlots.jl/stable/man/getting_started/)
+- in [PhyloNetworks](http://juliaphylo.github.io/PhyloNetworks.jl/latest/)'s documentation, see
+  * basic examples
+    [here](https://juliaphylo.github.io/PhyloNetworks.jl/dev/man/net_plot/#Network-Visualization)
+  * examples to annotate edges
+    [here](https://juliaphylo.github.io/PhyloNetworks.jl/dev/man/network_support/#support-for-tree-edges)
+- visualization with data in this
+  [tutorial](https://cecileane.github.io/networkPCM-workshop/topic6-visualization.html)
diff --git a/appveyor.yml b/appveyor.yml
deleted file mode 100644
index 2986315..0000000
--- a/appveyor.yml
+++ /dev/null
@@ -1,29 +0,0 @@
-environment:
-  matrix:
-  - julia_version: 1
-
-platform:
-  - x86 # 32-bit
-  - x64 # 64-bit
-
-branches:
-  only:
-    - master
-    - /release-.*/
-
-notifications:
-  - provider: Email
-    on_build_success: false
-    on_build_failure: false
-    on_build_status_changed: false
-
-install:
-  - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1"))
-
-build_script:
-  - echo "%JL_BUILD_SCRIPT%"
-  - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%"
-
-test_script:
-  - echo "%JL_TEST_SCRIPT%"
-  - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%"
diff --git a/docs/Project.toml b/docs/Project.toml
index 3c017cc..62d8675 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 PhyloNetworks = "33ad39ac-ed31-50eb-9b15-43d0656eaa72"
 PhyloPlots = "c0d5b6db-e3fc-52bc-a87d-1d050989ed3b"
 RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
diff --git a/docs/make.jl b/docs/make.jl
index f8fc1d8..2a965db 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,9 +3,16 @@ using Documenter
 using Pkg
 Pkg.add(PackageSpec(name="PhyloNetworks", rev="master"))
 
+using DocumenterInterLinks
+links = InterLinks(
+    "PhyloNetworks" => "https://juliaphylo.github.io/PhyloNetworks.jl/stable/objects.inv"
+)
 using PhyloNetworks
 using PhyloPlots
-DocMeta.setdocmeta!(PhyloPlots, :DocTestSetup, :(using PhyloPlots); recursive=true)
+# default loading of interlinked packages in all docstring examples
+DocMeta.setdocmeta!(PhyloPlots, :DocTestSetup,
+    :(using PhyloNetworks, PhyloPlots);
+    recursive=true)
 
 makedocs(
     sitename = "PhyloPlots.jl",
diff --git a/docs/src/man/adding_data.md b/docs/src/man/adding_data.md
index ec13e77..e141ca9 100644
--- a/docs/src/man/adding_data.md
+++ b/docs/src/man/adding_data.md
@@ -20,7 +20,7 @@ To add labels on edges (or nodes), we need to know their numbers. We can use the
 ```@example adding_data
 R"svg"(figname("adding_data1.svg"), width=6, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]); R"layout"([1 2]); # hide
-net = readTopology("(A,((B,#H1),((C)#H1, D)));") # hide
+net = readnewick("(A,((B,#H1),((C)#H1, D)));") # hide
 plot(net, showedgenumber=true);
 plot(net, showedgenumber=true, edgenumbercolor="red4");
 R"dev.off()" # hide
@@ -40,7 +40,7 @@ node), and the label that goes on it, like this:
 | 2      | "edge # 2" |
 
 After including the DataFrames package, we can define it as so:
-```@repl
+```@repl adding_data
 using DataFrames
 DataFrame(number=[1,2], label=["edge number 1","edge # 2"])
 ```
@@ -49,7 +49,7 @@ puts the text on the correct edges:
 ```@example adding_data
 R"svg"(figname("edge_labels_example.svg"), width=4, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
-net = readTopology("(A,((B,#H1),(C,(D)#H1)));") # hide
+net = readnewick("(A,((B,#H1),(C,(D)#H1)));") # hide
 plot(net, edgelabel=DataFrame(number = [1,2],
                               label = ["edge number 1", "edge # 2"]),
           edgelabelcolor="orangered", edgecex=[0.9,1.1]);
@@ -63,12 +63,12 @@ nothing # hide
 We can use the return values of [`plot`](@ref) to get information on the coordinates of
 different elements of the plot. Using this, we can add any other information we want.
 
-The [`plot`](@ref) function returns the following tuple:
+The [`plot`](@ref) function returns the following named tuple:
 ```
-(xmin, xmax, ymin, ymax,
- node_x, node_y, node_yB, node_yE,
- edge_xB, edge_xE, edge_yB, edge_yE,
- nodedataframe, edgedataframe)
+(:xmin, :xmax, :ymin, :ymax,
+ :node_x,    :node_y,    :node_y_lo, :node_y_hi,
+ :edge_x_lo, :edge_x_hi, :edge_y_lo, :edge_y_hi,
+ :node_data, :edge_data)
 ```
 See the documentation for descriptions of these elements: [`plot`](@ref)
 
@@ -79,7 +79,7 @@ Here's example code that adds bars to denote clades in the margin:
 ```@example adding_data
 R"svg"(figname("side_bars.svg"), width=4, height=4) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
-net = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
+net = readnewick("(((((((t1,t2),t3),t4),t5),(t6,t7)),(t8,t9)),t10);");
 plot(net, xlim=(1,10))
 using RCall # to send any R command, to make further plot modifications
 R"segments"([9, 9, 9], [0.8, 7.8, 9.8], [9, 9, 9], [7.2, 9.2, 10.2])
@@ -97,7 +97,7 @@ because they contain the default range of the x axis; `xmin` and `xmax`.
 
 ```@example adding_data
 res = plot(net);
-res[1:2]
+res[[:xmin,:xmax]]
 ```
 
 Looking at `xmin` and `xmax` returned by default, we can see that the x
@@ -120,3 +120,17 @@ using RCall # add (install) the RCall package prior to 'using' it
 R"segments"([9, 9, 9], [0.8, 7.8, 9.8], [9, 9, 9], [7.2, 9.2, 10.2])
 R"text"([9.5, 9.5, 9.5], [4, 8.5, 10], ["C", "B", "A"])
 ```
+
+# Beyond
+
+To go beyond, we can access data on the node & edges to use them as we wish.
+We can access the coordinates of points & segments and more data like this:
+
+```@repl adding_data
+res[:node_x] # x coordinate. similarly try res[:node_y]
+hcat(res[:node_y_lo], res[:node_y_hi])
+DataFrame(edge_x_lo=res[:edge_x_lo], edge_x_hi=res[:edge_x_hi],
+          edge_y_lo=res[:edge_y_lo], edge_y_hi=res[:edge_y_hi])
+res[:node_data]
+res[:edge_data]
+```
diff --git a/docs/src/man/better_edges.md b/docs/src/man/better_edges.md
index 2333bcc..580147d 100644
--- a/docs/src/man/better_edges.md
+++ b/docs/src/man/better_edges.md
@@ -15,7 +15,7 @@ but by switching it to `:majortree`, we can draw minor hybrid edges as diagonal
 ```@example better_edges
 R"svg"(figname("style_example.svg"), width=3, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
-net = readTopology("(A,((B,#H1),(C,(D)#H1)));") # hide
+net = readnewick("(A,((B,#H1),(C,(D)#H1)));") # hide
 plot(net, style=:majortree);
 R"dev.off()" # hide
 nothing # hide
@@ -31,7 +31,7 @@ lines. For this, we'll use a network that has branch lengths:
 R"svg"(figname("edge_len_example.svg"), width=6, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
 R"layout"([1 2]) # hide
-net = readTopology("(A:3.3,((B:1.5,#H1:0.5):1.5,((C:1)#H1:1.8,D:1.1):.2):0.3);")
+net = readnewick("(A:3.3,((B:1.5,#H1:0.5):1.5,((C:1)#H1:1.8,D:1.1):.2):0.3);")
 df = DataFrame(number=[-3,3], label=["N","H1"]); # hide
 plot(net, useedgelength=true, ylim = [-1, 5.5], nodelabel = df); # hide
 R"text"([3], [0], ["useedgelength=true"]) # hide
@@ -73,8 +73,8 @@ Time-inconsistent networks like these ones below might cause confusion:
 R"svg"(figname("edge_len_example2.svg"), width=6, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
 R"layout"([1 2]) # hide
-net1 = readTopology("(A:3.3,((B:1.5,#H1:1.2):1.5,((C:1.8)#H1:1,D:1.1):.2):0.3);");
-net2 = readTopology("(A:3.3,((B:1.5,#H1:0.2):1.5,((C:1)#H1:1.8,D:1.1):.2):0.3);");
+net1 = readnewick("(A:3.3,((B:1.5,#H1:1.2):1.5,((C:1.8)#H1:1,D:1.1):.2):0.3);");
+net2 = readnewick("(A:3.3,((B:1.5,#H1:0.2):1.5,((C:1)#H1:1.8,D:1.1):.2):0.3);");
 plot(net1, useedgelength=true); # hide
 plot(net2, useedgelength=true); # hide
 R"dev.off()" # hide
diff --git a/docs/src/man/getting_started.md b/docs/src/man/getting_started.md
index 926770c..c3f50c1 100644
--- a/docs/src/man/getting_started.md
+++ b/docs/src/man/getting_started.md
@@ -16,14 +16,14 @@ using PhyloPlots
 ```
 Then read the topology
 ```@repl getting_started
-net = readTopology("(A,((B,#H1),(C,(D)#H1)));")
+net = readnewick("(A,((B,#H1),(C,(D)#H1)));")
 ```
 and call `plot`, as shown below.
 
 ```@example getting_started
 R"svg"(figname("gettingstarted.svg"), width=3, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
-net = readTopology("(A,((B,#H1),(C,(D)#H1)));") # hide
+net = readnewick("(A,((B,#H1),(C,(D)#H1)));") # hide
 plot(net);
 R"dev.off()" # hide
 nothing # hide
diff --git a/docs/src/man/untangling_edges.md b/docs/src/man/untangling_edges.md
index 85be8f3..d33b3db 100644
--- a/docs/src/man/untangling_edges.md
+++ b/docs/src/man/untangling_edges.md
@@ -20,7 +20,7 @@ whose child edges we should rotate.
 ```@example untangling
 R"svg"(figname("untangling1.svg"), width=3, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
-net = readTopology("(A,((B,#H1),(C,(D)#H1)));") # hide
+net = readnewick("(A,((B,#H1),(C,(D)#H1)));") # hide
 plot(net, shownodenumber=true);
 R"dev.off()" # hide
 nothing # hide
@@ -32,7 +32,7 @@ As we can see, rotating edges around node `-5` will make for a prettier network.
 ```@example untangling
 R"svg"(figname("untangling2.svg"), width=3, height=3) # hide
 R"par"(mar=[.1,.1,.1,.1]) # hide
-net = readTopology("(A,((B,#H1),(C,(D)#H1)));") # hide
+net = readnewick("(A,((B,#H1),(C,(D)#H1)));") # hide
 rotate!(net, -5)
 plot(net)
 R"dev.off()" # hide
diff --git a/src/PhyloPlots.jl b/src/PhyloPlots.jl
index f1d5c3a..98fe265 100644
--- a/src/PhyloPlots.jl
+++ b/src/PhyloPlots.jl
@@ -18,7 +18,6 @@ export sexp
 
 include("phylonetworksPlots.jl")
 include("plotRCall.jl")
-include("substitutionmodels.jl")
 include("rexport.jl")
 
 end # of module
diff --git a/src/phylonetworksPlots.jl b/src/phylonetworksPlots.jl
index 9fc70e3..85dd586 100644
--- a/src/phylonetworksPlots.jl
+++ b/src/phylonetworksPlots.jl
@@ -7,7 +7,7 @@
 Calculate coordinates for plotting later with Gadfly or RCall.
 
 Actually modifies some (minor) attributes of the network,
-as it calls `directEdges!` and `preorder!`.
+as it calls `directedges!` and `preorder!`.
 
 output: tuple with the following elements, in which the order of
 nodes corresponds to the order in `net.node`, and the order of
@@ -38,39 +38,43 @@ edges corresponds to that in `net.edge` (filtered to minor edges as needed).
    the same as that of the mid-point, given by `node_x`.
 9. `minoredge_xB`: x coordinate for the Beginning and ...
 10. `minoredge_xE`: ... End of the diagonal segment of each minor hybrid edge,
-   in the same order as in `filter(e -> !e.isMajor, net.edge)`.
+   in the same order as in `filter(e -> !e.ismajor, net.edge)`.
 11. `minoredge_yB`: y coordinate for the beginning and ...
 12. `minoredge_yE`: ... end of the diagonal segment of each minor hybrid edge.
 13-16. `xmin`, `xmax`, `ymin`, `ymax`: ranges for the x and y axes.
 """
-function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimpleHybridLines::Bool)
+function edgenode_coordinates(
+    net::HybridNetwork,
+    useedgelength::Bool,
+    useSimpleHybridLines::Bool
+)
     try
-        directEdges!(net)   # to update isChild1
+        directedges!(net)   # to update ischild1
     catch e
         if isa(e, PhyloNetworks.RootMismatch)
             e = PhyloNetworks.RootMismatch( e.msg * "\nPlease change the root, perhaps using rootatnode! or rootatedge!")
         end
         rethrow(e)
     end
-    preorder!(net)       # to update net.nodes_changed: true pre-ordering
+    preorder!(net)       # to update net.vec_node: true pre-ordering
 
     # determine y for each node = y of its parent edge: post-order traversal
     # also [yB,yE] for each internal node: range of y's of all children nodes
-    # y max is the numTaxa + number of minor edges
+    # y max is the numtaxa + number of minor edges
     ymin = 1.0;
-    ymax = net.numTaxa
+    ymax = net.numtaxa
     if !useSimpleHybridLines
-        ymax += sum(!e.isMajor for e in net.edge)
+        ymax += sum(!e.ismajor for e in net.edge)
     end
 
-    node_y  = zeros(Float64, net.numNodes) # order: in net.nodes, *!not in nodes_changed!*
-    node_yB = zeros(Float64,net.numNodes) # min (B=begin) and max (E=end)
-    node_yE = zeros(Float64,net.numNodes) #   of at children's nodes
-    edge_yB = zeros(Float64,net.numEdges) # yE of edge = y of child node
+    node_y  = zeros(Float64, net.numnodes) # order: in net.nodes, *!not in vec_node!*
+    node_yB = zeros(Float64,net.numnodes) # min (B=begin) and max (E=end)
+    node_yE = zeros(Float64,net.numnodes) #   of at children's nodes
+    edge_yB = zeros(Float64,net.numedges) # yE of edge = y of child node
     # set node_y of leaves: follow cladewise order
     # also sets edge_yB of minor hybrid edges
     nexty = ymax # first tips at the top, last at bottom
-    cladewise_queue = copy(net.node[net.root].edge) # the child edges of root
+    cladewise_queue = copy(getroot(net).edge) # the child edges of root
     # print("queued the root's children's indices: "); @show queue
     while !isempty(cladewise_queue)
         cur_edge = pop!(cladewise_queue); # deliberate choice over shift! for cladewise order
@@ -86,13 +90,13 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
 
         # only for new hybrid lines:
         # increment spacing and add to edge_yB if parent edge is minor
-        if !cur_edge.isMajor && !useSimpleHybridLines
+        if !cur_edge.ismajor && !useSimpleHybridLines
             edge_yB[findfirst(x->x===cur_edge, net.edge)] = nexty
             nexty -= 1
         end
 
         # push children edges if this is a major edge:
-        if cur_edge.isMajor
+        if cur_edge.ismajor
             for e in cur_child.edge
                 if getparent(e) === cur_child # don't go backwards
                     push!(cladewise_queue, e)
@@ -103,7 +107,7 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
 
     # set node_y of internal nodes: follow post-order
     for i=length(net.node):-1:1
-        nn = net.nodes_changed[i]
+        nn = net.vec_node[i]
         !nn.leaf || continue # previous loop took care of leaves
         ni = findfirst(x -> x===nn, net.node)
         node_yB[ni]=ymax; node_yE[ni]=ymin;
@@ -113,12 +117,12 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
             if nn == getparent(e) # if e = child of node
                 if useSimpleHybridLines
                     # old simple hybrid lines
-                    if e.isMajor || nomajorchild
+                    if e.ismajor || nomajorchild
                         cc = getchild(e)
                         yy = node_y[findfirst(x -> x===cc, net.node)]
                         yy!==nothing || error("oops, child $(cc.number) has not been visited before node $(nn.number).")
                     end
-                    if e.isMajor
+                    if e.ismajor
                         nomajorchild = false # we found a child edge that is a major edge
                         node_yB[ni] = min(node_yB[ni], yy)
                         node_yE[ni] = max(node_yE[ni], yy)
@@ -128,7 +132,7 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
                     end
                 else
                     # new pretty hybrid lines
-                    if e.isMajor
+                    if e.ismajor
                         cc = getchild(e)
                         child_y = node_y[findfirst(x -> x===cc, net.node)]
                         child_y!==nothing || error("oops, child $(cc.number) has not been visited before node $(nn.number).")
@@ -176,13 +180,13 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
         # setting elen such that the age of each node = 1 + age of oldest child
         # (including minor hybrid edges): need true post-ordering.
         # calculating node ages first, elen will be calculated later.
-        elen     = zeros(Float64,net.numEdges)
-        node_age = zeros(Float64,net.numNodes)
+        elen     = zeros(Float64,net.numedges)
+        node_age = zeros(Float64,net.numnodes)
         for i=length(net.node):-1:1 # post-order traversal
-            if net.nodes_changed[i].leaf continue; end
-            ni = findfirst(x -> x===net.nodes_changed[i], net.node)
-            for e in net.nodes_changed[i].edge # loop over children only
-                if net.nodes_changed[i] == (e.isChild1 ? e.node[2] : e.node[1])
+            if net.vec_node[i].leaf continue; end
+            ni = findfirst(x -> x===net.vec_node[i], net.node)
+            for e in net.vec_node[i].edge # loop over children only
+                if net.vec_node[i] == (e.ischild1 ? e.node[2] : e.node[1])
                     node_age[ni] = max(node_age[ni], 1 +
                      node_age[findfirst(x -> x=== getchild(e), net.node)])
                 end
@@ -197,15 +201,15 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
     # determine xB,xE for each edge: pre-order traversal, uses branch lengths
     # then x and yB,yE for each node: x=xE of parent edge
     xmin = 1.0; xmax=xmin
-    node_x  = zeros(Float64,net.numNodes) # order: in net.nodes, *!not in nodes_changed!*
-    edge_xB = zeros(Float64,net.numEdges) # min (B=begin) and max (E=end)
-    edge_xE = zeros(Float64,net.numEdges) # xE-xB = edge length
-    node_x[net.root] = xmin # root node: x=xmin=0
+    node_x  = zeros(Float64,net.numnodes) # order: in net.nodes, *!not in vec_node!*
+    edge_xB = zeros(Float64,net.numedges) # min (B=begin) and max (E=end)
+    edge_xE = zeros(Float64,net.numedges) # xE-xB = edge length
+    node_x[net.rooti] = xmin # root node: x=xmin=0
     for i=2:length(net.node)              # true pre-order, skipping the root (i=1)
-        ni = findfirst(x -> x===net.nodes_changed[i], net.node)
+        ni = findfirst(x -> x===net.vec_node[i], net.node)
         ei = nothing # index of major parent edge of current node
-        for e in net.nodes_changed[i].edge
-            if e.isMajor && net.nodes_changed[i] == e.node[e.isChild1 ? 1 : 2] # major parent edge
+        for e in net.vec_node[i].edge
+            if e.ismajor && net.vec_node[i] == e.node[e.ischild1 ? 1 : 2] # major parent edge
                 ei = findfirst(x -> x===e, net.edge)
                 break
             end
@@ -228,8 +232,8 @@ function edgenode_coordinates(net::HybridNetwork, useedgelength::Bool, useSimple
     minoredge_yB = Float64[]
     minoredge_yE = Float64[]
 
-    for i=1:net.numEdges
-        if !net.edge[i].isMajor # minor hybrid edges
+    for i=1:net.numedges
+        if !net.edge[i].ismajor # minor hybrid edges
             # indices of child and parent nodes
             cni = findfirst(x -> x===getchild( net.edge[i]), net.node)
             pni = findfirst(x -> x===getparent(net.edge[i]), net.node)
@@ -270,7 +274,10 @@ Check data frame for node annotations:
 - remove rows with no node numbers
 - warning if some node numbers in the data are not in the network.
 """
-function check_nodedataframe(net::HybridNetwork, nodelabel::DataFrame)
+function check_nodedataframe(
+    net::HybridNetwork,
+    nodelabel::DataFrame
+)
     labelnodes = size(nodelabel,1)>0
     if (labelnodes && (size(nodelabel,2)<2 ||
             !(nonmissingtype(eltype(nodelabel[!,1])) <: Integer)))
@@ -309,16 +316,22 @@ Columns of output data frame:
 - lab: node label
 - lea: is leaf?
 """
-function prepare_nodedataframe(net::HybridNetwork, nodelabel::DataFrame,
-        shownodenumber::Bool, shownodename::Bool, labelnodes::Bool,
-        node_x::Array{Float64,1}, node_y::Array{Float64,1})
-    nrows = (shownodenumber || shownodename || labelnodes ? net.numNodes : net.numTaxa)
+function prepare_nodedataframe(
+    net::HybridNetwork,
+    nodelabel::AbstractDataFrame,
+    shownodenumber::Bool,
+    shownodename::Bool,
+    labelnodes::Bool,
+    node_x::Array{Float64,1},
+    node_y::Array{Float64,1}
+)
+    nrows = (shownodenumber || shownodename || labelnodes ? net.numnodes : net.numtaxa)
     ndf = DataFrame(:name => Vector{String}(undef,nrows),
         :num => Vector{String}(undef,nrows), :lab => Vector{String}(undef,nrows),
         :lea => Vector{Bool}(  undef,nrows), :x => Vector{Float64}( undef,nrows),
         :y => Vector{Float64}( undef,nrows), copycols=false)
     j=1
-    for i=1:net.numNodes
+    for i=1:net.numnodes
     if net.node[i].leaf  || shownodenumber || shownodename || labelnodes
         ndf[j,:name] = net.node[i].name
         ndf[j,:num] = string(net.node[i].number)
@@ -354,12 +367,20 @@ Return data frame with columns
 - hyb: is hybrid?
 - min: is minor?
 """
-function prepare_edgedataframe(net::HybridNetwork, edgelabel::DataFrame, style::Symbol,
-        edge_xB::Array{Float64,1}, edge_xE::Array{Float64,1},
-        edge_yB::Array{Float64,1}, edge_yE::Array{Float64,1},
-        minoredge_xB::Array{Float64,1}, minoredge_xE::Array{Float64,1},
-        minoredge_yB::Array{Float64,1}, minoredge_yE::Array{Float64,1})
-    nrows = net.numEdges
+function prepare_edgedataframe(
+    net::HybridNetwork,
+    edgelabel::AbstractDataFrame,
+    style::Symbol,
+    edge_xB::Array{Float64,1},
+    edge_xE::Array{Float64,1},
+    edge_yB::Array{Float64,1},
+    edge_yE::Array{Float64,1},
+    minoredge_xB::Array{Float64,1},
+    minoredge_xE::Array{Float64,1},
+    minoredge_yB::Array{Float64,1},
+    minoredge_yE::Array{Float64,1}
+)
+    nrows = net.numedges
     edf = DataFrame(:len => Vector{String}(undef,nrows),
         :gam => Vector{String}(undef,nrows), :num => Vector{String}(undef,nrows),
         :lab => Vector{String}(undef,nrows), :hyb => Vector{Bool}(undef,nrows),
@@ -384,7 +405,7 @@ function prepare_edgedataframe(net::HybridNetwork, edgelabel::DataFrame, style::
       end
     end
     j=1   # index of row in edf
-    imh=1 # index of minor hybrid edge in filter(ee -> !ee.isMajor, net.edge)
+    imh=1 # index of minor hybrid edge in filter(ee -> !ee.ismajor, net.edge)
     for i = 1:length(net.edge)
         ee = net.edge[i]
         edf[j,:len] = (ee.length==-1.0 ? "" : @sprintf("%0.3g",ee.length))
@@ -398,10 +419,10 @@ function prepare_edgedataframe(net::HybridNetwork, edgelabel::DataFrame, style::
                 @sprintf("%0.3g",edgelabel[je,2]) : string(edgelabel[je,2])))
         end
         edf[j,:hyb] = ee.hybrid
-        edf[j,:min] = !ee.isMajor
+        edf[j,:min] = !ee.ismajor
         edf[j,:y] = (edge_yB[i] + edge_yE[i])/2
         edf[j,:x] = (edge_xB[i] + edge_xE[i])/2
-        if style == :majortree && !ee.isMajor
+        if style == :majortree && !ee.ismajor
             edf[j,:y] = (minoredge_yB[imh] + minoredge_yE[imh])/2
             edf[j,:x] = (minoredge_xB[imh] + minoredge_xE[imh])/2
             imh += 1
diff --git a/src/plotRCall.jl b/src/plotRCall.jl
index 702c9d6..6faee2c 100755
--- a/src/plotRCall.jl
+++ b/src/plotRCall.jl
@@ -1,135 +1,123 @@
-@deprecate plot(net::HybridNetwork, ::Symbol; useEdgeLength=false::Bool,
-  mainTree=false::Bool, showTipLabel=true::Bool, showNodeNumber=false::Bool,
-  showEdgeLength=false::Bool, showGamma=false::Bool,
-  edgeColor="black"::String,
-  majorHybridEdgeColor="deepskyblue4"::String,
-  minorHybridEdgeColor="deepskyblue"::String,
-  showEdgeNumber=false::Bool, showIntNodeLabel=false::Bool,
-  edgeLabel=DataFrame()::DataFrame, nodeLabel=DataFrame()::DataFrame,
-  xlim=Float64[]::Array{Float64,1}, ylim=Float64[]::Array{Float64,1},
-  tipOffset=0.0::Float64, tipcex=1.0::Float64,
-  style=:fulltree::Symbol, arrowlen=(style==:majortree ? 0 : 0.1)::Real
-) plot(net::HybridNetwork; useedgelength=useEdgeLength,
-  showtiplabel=showTipLabel, shownodenumber=showNodeNumber,
-  showedgelength=showEdgeLength, showgamma=showGamma,
-  edgecolor=edgeColor,
-  majorhybridedgecolor=majorHybridEdgeColor,
-  minorhybridedgecolor=minorHybridEdgeColor,
-  showedgenumber=showEdgeNumber, shownodelabel=showIntNodeLabel,
-  edgelabel=edgeLabel, nodelabel=nodeLabel,
-  xlim=xlim, ylim=ylim,
-  tipoffset=tipOffset, tipcex=tipcex,
-  style=style, arrowlen=arrowlen)
-
 """
     plot(net::HybridNetwork)
 
-Plot a network using R graphics. Optional arguments are listed below.
+Plot a network with edges going from left to right and taxa (leaves) placed on
+the right, using R graphics. Optional arguments are listed below.
 
 ## lines forming the network:
 
-- `useedgelength = false` : if true, the tree edges and major hybrid edges are
+- `useedgelength = false`: if true, the tree edges and major hybrid edges are
   drawn proportionally to their length. Minor hybrid edges are not, however.
   Note that edge lengths in coalescent units may scale very poorly with time.
-- `style = :fulltree` : symbol indicating the style of the diagram
+- `style = :fulltree`: symbol indicating the style of the diagram
   * `:majortree` will simply draw minor edges onto the major tree.
   * `:fulltree` will draw minor edges as their own branches in the tree,
     in the same style used by [icytree](https://icytree.org). This is
     useful for overlapping or confusing networks.
-- `arrowlen` : the length of the arrow tips in the full tree style. if `style = :fulltree`, then
-  `arrowlen = 0.2`. otherwise, `arrowlen = 0`, which makes the arrows appear as segments.
+- `arrowlen`: the length of the arrow tips in the full tree style.
+  The default is 0.1 if `style = :fulltree`,
+  and 0 if `style = :majortree` (making the arrows appear as segments).
 - `edgewidth=1`: width of horizontal (not diagonal) edges. To vary them,
-  use a dictionary to map the number of each edge to its desired witdth.
-- `xlim`, `ylim` : array of 2 values, to determine the axes limits
+  use a dictionary to map the number of each edge to its desired width.
+- `xlim`, `ylim`: array of 2 values, to determine the axes limits.
 
 ## tip annotations:
 
-- `showtiplabel = true` : if true, taxon labels (names) are shown.
-- `tipoffset = 0.0` : to offset tip labels
-- `tipcex=1.0`: character expansion for tip and internal node names
+- `showtiplabel = true`: if true, taxon labels (names) are shown.
+- `tipoffset = 0`: to offset tip labels.
+- `tipcex = 1`: character expansion for tip and internal node names.
 
 ## nodes & edges annotations:
 
-- `shownodelabel = false` : if true, internal nodes are labelled with their names.
+- `shownodelabel = false`: if true, internal nodes are labelled with their names.
   Useful for hybrid nodes, which do have tags like 'H1'.
-- `shownodenumber = false` : if true, nodes are labelled with the number used internally.
-- `showedgenumber = false` : if true, edges are labelled with the number used internally.
-- `showedgelength = false` : if true, edges are labelled with their length (above)
-- `showgamma = false` : if true, hybrid edges are labelled with their heritability (below)
-- `edgelabel = DataFrame()` : dataframe with two columns: the first with edge numbers,
+- `shownodenumber = false`: if true, nodes are labelled with the number used internally.
+- `showedgenumber = false`: if true, edges are labelled with the number used internally.
+- `showedgelength = false`: if true, edges are labelled with their length (above).
+- `showgamma = false`: if true, hybrid edges are labelled with their heritability (below).
+- `edgelabel = DataFrame()`: dataframe with two columns: the first with edge numbers,
   the second with labels (like bootstrap values) to annotate edges. empty by default.
-- `nodelabel = DataFrame()` : dataframe with two columns: the first with node numbers,
+- `nodelabel = DataFrame()`: dataframe with two columns: the first with node numbers,
   the second with labels (like bootstrap values for hybrid relationships)
   to annotate nodes. empty by default.
-- `nodecex=1.0`: character expansion for labels in the `nodelabel` data frame
-- `edgecex=1.0`: character expansion for labels in the `edgelabel` data frame
+- `nodecex = 1`: character expansion for labels in the `nodelabel` data frame.
+- `edgecex = 1`: character expansion for labels in the `edgelabel` data frame.
 
 ## colors:
 
-- `edgecolor = "black"` : color for tree edges.
-- `majorhybridedgecolor = "deepskyblue4"` : color for major hybrid edges
-- `minorhybridedgecolor = "deepskyblue"` : color for minor hybrid edges
-- `edgenumbercolor = "grey"`: color for edge numbers
-- `edgelabelcolor = "black"`: color for labels in the `edgelabel` data frame
-- `nodelabelcolor = "black"`: color for labels in the `nodelabel` data frame
+- `edgecolor = "black"`: color for tree edges.
+- `majorhybridedgecolor = "deepskyblue4"`: color for major hybrid edges.
+- `minorhybridedgecolor = "deepskyblue"`: color for minor hybrid edges.
+- `edgenumbercolor = "grey"`: color for edge numbers.
+- `edgelabelcolor = "black"`: color for labels in the `edgelabel` data frame.
+- `nodelabelcolor = "black"`: color for labels in the `nodelabel` data frame.
 
-Output the following tuple, that can be used for downstream plot annotations
+Output the following named tuple, that can be used for downstream plot annotations
 with RCall:
 
 ```
 (xmin, xmax, ymin, ymax,
- node_x, node_y, node_yB, node_yE,
- edge_xB, edge_xE, edge_yB, edge_yE,
- ndf, edf)
+ node_x,    node_y,    node_y_lo, node_y_hi,
+ edge_x_lo, edge_x_hi, edge_y_lo, edge_y_hi,
+ node_data, edge_data)
 ```
 
-1. `xmin` : minimum x value of the plot
-2. `xmax` : maximum x value of the plot
-3. `ymin` : minimum y value of the plot
-4. `ymax` : maximum y value of the plot
-5. `node_x` : x values of the nodes in net.node in their respective order
-6. `node_y` : y values of the nodes
-7. `node_yB` : y value of the beginning of the verticle bar
-8. `node_yE` : y value of the end of the verticle bar
-9. `edge_xB` : x value of the beginning of the edges in net.edge in their respective order
-10. `edge_xE` : x value of the end of the edges
-11. `edge_yB` : y value of the beginning of the edges
-12. `edge_yE` : y value of the end of the edges
-13. `ndf` : node data frame: see section [Adding labels](@ref) for more
-14. `edf` : edge data frame
+1. `:xmin`: minimum x value of the plot
+2. `:xmax`: maximum x value of the plot
+3. `:ymin`: minimum y value of the plot
+4. `:ymax`: maximum y value of the plot
+5. `:node_x`: x values of the nodes in net.node in their respective order
+6. `:node_y`: y values of the nodes
+7. `:node_y_lo`: y value of the beginning of the vertical bar representing the clade at each node
+8. `:node_y_hi`: y value of the end of the vertical bar
+9. `:edge_x_lo`: x value of the beginning of the edges in `net.edge` in their respective order
+10. `:edge_x_hi`: x value of the end of the edges
+11. `:edge_y_lo`: y value of the beginning of the edges
+12. `:edge_y_hi`: y value of the end of the edges
+13. `:node_data`: node data frame: see section [Adding labels](@ref) for more
+14. `:edge_data`: edge data frame
 
 Note that `plot` actually modifies some (minor) attributes of the network,
-as it calls `directEdges!` and `preorder!`.
+as it calls `PhyloNetworks.directedges!` and `PhyloNetworks.preorder!`.
 
 If hybrid edges cross tree and major edges, you may choose to rotate some tree
-edges to eliminate crossing edges, using `rotate!`
-(in [`PhyloNetworks`](http://juliaphylo.github.io/PhyloNetworks.jl/latest/lib/public/#PhyloNetworks.rotate!)).
+edges to eliminate crossing edges, using
+[`PhyloNetworks.rotate!`](https://juliaphylo.github.io/PhyloNetworks.jl/dev/lib/public/#PhyloNetworks.rotate!-Tuple%7BHybridNetwork,%20Integer%7D).
 
 **Alternative**: a tree or network can be exported with [`sexp`](@ref)
 and then displayed with R's "plot" and all its options.
 """
-function plot(net::HybridNetwork; useedgelength=false::Bool,
-    showtiplabel=true::Bool, shownodenumber=false::Bool,
-    showedgelength=false::Bool, showgamma=false::Bool,
-    edgecolor="black"::String,
-    majorhybridedgecolor="deepskyblue4"::String,
-    minorhybridedgecolor="deepskyblue"::String,
-    showedgenumber=false::Bool, shownodelabel=false::Bool,
-    edgelabel=DataFrame()::DataFrame, nodelabel=DataFrame()::DataFrame,
-    xlim=Float64[]::Array{Float64,1}, ylim=Float64[]::Array{Float64,1},
-    tipoffset=0.0::Float64, tipcex=1.0::Float64,
-    nodecex=1.0::Float64, edgecex=1.0::Float64,
-    style=:fulltree::Symbol, arrowlen=(style==:majortree ? 0 : 0.1)::Real,
+function plot(
+    net::HybridNetwork;
+    useedgelength::Bool=false,
+    showtiplabel::Bool=true,
+    shownodenumber::Bool=false,
+    showedgelength::Bool=false,
+    showgamma::Bool=false,
+    edgecolor::AbstractString="black",
+    majorhybridedgecolor::AbstractString="deepskyblue4",
+    minorhybridedgecolor::AbstractString="deepskyblue",
+    showedgenumber::Bool=false,
+    shownodelabel::Bool=false,
+    edgelabel::AbstractDataFrame=DataFrame(),
+    nodelabel::AbstractDataFrame=DataFrame(),
+    xlim = nothing,
+    ylim = nothing,
+    tipoffset = 0,
+    tipcex = 1,
+    nodecex = 1,
+    edgecex = 1,
+    style::Symbol=:fulltree,
+    arrowlen::Real=(style==:majortree ? 0 : 0.1),
     edgewidth = 1,
     edgenumbercolor = "grey", # don't limit the type because R accepts many types
     edgelabelcolor = "black", # and these colors are used as is
     nodelabelcolor = "black",
 )
-
-    if net.node[net.root].leaf
+    if getroot(net).leaf
         @warn """The network is rooted at a leaf: the plot won't look good.
             Try rooting the network on the edge adjacent to that leaf, with
-            rootonedge!(network_name, $(net.node[net.root].edge[1].number))"""
+            rootonedge!(network_name, $(getroot(net).edge[1].number))"""
     end
     (edge_xB, edge_xE, edge_yB, edge_yE, node_x, node_y, node_yB, node_yE,
      hybridedge_xB, hybridedge_xE, hybridedge_yB, hybridedge_yE,
@@ -146,16 +134,20 @@ function plot(net::HybridNetwork; useedgelength=false::Bool,
         ymax += expfacy
     end
     xmax += tipoffset
-    if length(xlim)==2
+    if !isnothing(xlim)
+        length(xlim) == 2 ||
+          error("xlim needs to contain 2 values: lower and upper limits. defaults: [$xmin,$xmax]")
         xmin=xlim[1]; xmax=xlim[2]
     end
-    if length(ylim)==2
+    if !isnothing(ylim)
+        length(ylim) == 2 ||
+          error("ylim needs to contain 2 values: lower and upper limits. defaults: [$ymin,$ymax]")
         ymin=ylim[1]; ymax=ylim[2]
     end
     leaves = [n.leaf for n in net.node]
     eCol = fill(edgecolor, length(net.edge))
     eCol[ [ e.hybrid  for e in net.edge] ] .= majorhybridedgecolor
-    eCol[ [!e.isMajor for e in net.edge] ] .= minorhybridedgecolor
+    eCol[ [!e.ismajor for e in net.edge] ] .= minorhybridedgecolor
 
     if isa(edgewidth, Number)
       edgewidth_vec = edgewidth
@@ -187,7 +179,7 @@ function plot(net::HybridNetwork; useedgelength=false::Bool,
     R"segments"(node_x, node_yB, node_x, node_yE, col=edgecolor)
     if showtiplabel
       R"text"(node_x[leaves] .+ tipoffset, node_y[leaves],
-              tipLabels(net), adj=0, font=3, cex=tipcex)
+              tiplabels(net), adj=0, font=3, cex=tipcex)
     end
     if shownodelabel
       R"text"(ndf[.!ndf[!,:lea],:x], ndf[.!ndf[!,:lea],:y],
@@ -210,7 +202,7 @@ function plot(net::HybridNetwork; useedgelength=false::Bool,
     if showedgelength
       R"text"(edf[!,:x], edf[!,:y], edf[!,:len], adj=[.5,1.])
     end
-    if showgamma && net.numHybrids>0
+    if showgamma && net.numhybrids>0
       im = edf[!,:hyb] .& edf[!,:min]
       iM = edf[!,:hyb] .& .!edf[!,:min]
       R"text"(edf[im,:x], edf[im,:y], edf[im,:gam],
@@ -221,6 +213,10 @@ function plot(net::HybridNetwork; useedgelength=false::Bool,
     if showedgenumber
       R"text"(edf[!,:x], edf[!,:y], edf[!,:num], adj=[.5,0], col=edgenumbercolor)
     end
-    return (xmin, xmax, ymin, ymax, node_x, node_y, node_yB, node_yE,
-      edge_xB, edge_xE, edge_yB, edge_yE, ndf, edf)
+    return (xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
+      node_x=node_x, node_y=node_y,
+      node_y_lo=node_yB, node_y_hi=node_yE,
+      edge_x_lo=edge_xB, edge_x_hi=edge_xE,
+      edge_y_lo=edge_yB, edge_y_hi=edge_yE,
+      node_data=ndf, edge_data=edf)
 end
diff --git a/src/rexport.jl b/src/rexport.jl
index f8a7b6f..ee45b28 100644
--- a/src/rexport.jl
+++ b/src/rexport.jl
@@ -15,7 +15,7 @@ using plotting facilities in R.
 ```julia-repl
 julia> using RCall
 julia> using PhyloNetworks
-julia> net = readTopology("(((A:.2,(B:.1)#H1:.1::0.9):.1,(C:.11,#H1:.01::0.1):.19):.1,D:.4);");
+julia> net = readnewick("(((A:.2,(B:.1)#H1:.1::0.9):.1,(C:.11,#H1:.01::0.1):.19):.1,D:.4);");
 R> library(ape); # type $ to switch from julia to R
 R> $net
 
@@ -58,7 +58,7 @@ R> plot(net)
 ```
 """ #"
 function sexp(net::HybridNetwork)
-    PhyloNetworks.resetNodeNumbers!(net)
+    PhyloNetworks.resetnodenumbers!(net)
     ntips = length(net.leaf)
     totalnodes = length(net.node)
     Nnode = totalnodes - ntips
@@ -73,7 +73,7 @@ function sexp(net::HybridNetwork)
     if any(.!ismissing.(edgeLength))
         phy[Symbol("edge.length")] = edgeLength
     end
-    if net.numHybrids > 0
+    if net.numhybrids > 0
         reticulation = PhyloNetworks.minorreticulationmatrix(net)
         reticulationGamma = PhyloNetworks.minorreticulationgamma(net)
         reticulationLength = PhyloNetworks.minorreticulationlength(net)
@@ -86,7 +86,7 @@ function sexp(net::HybridNetwork)
         end
     end
     sobj = RCall.protect(sexp(phy)) # RObject
-    if net.numHybrids == 0
+    if net.numhybrids == 0
         setclass!(sobj, sexp("phylo"))
     else
         setclass!(sobj, sexp(["evonet", "phylo"]))
@@ -113,7 +113,7 @@ not exported: [`sexp`](@ref) is the best way to go.
 # Examples
 
 ```julia-repl
-julia> net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);");
+julia> net = readnewick("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);");
 
 julia> phy = PhyloPlots.rexport(net)
 RObject{VecSxp}
@@ -171,10 +171,10 @@ function rexport(net::HybridNetwork; maintree::Bool=false, useedgelength::Bool=t
 # worry about R object created within the function not accessible from outside:
 # can it be garbage collected?
 
-    if maintree && net.numHybrids > 0
-        net = majorTree(net)
+    if maintree && net.numhybrids > 0
+        net = majortree(net)
     end
-    PhyloNetworks.resetNodeNumbers!(net)
+    PhyloNetworks.resetnodenumbers!(net)
     ntips = length(net.leaf)
     totalnodes = length(net.node)
     Nnode = totalnodes - ntips
@@ -192,7 +192,7 @@ function rexport(net::HybridNetwork; maintree::Bool=false, useedgelength::Bool=t
             """
         end
     end
-    if net.numHybrids > 0
+    if net.numhybrids > 0
         reticulation = PhyloNetworks.minorreticulationmatrix(net)
         reticulationGamma = PhyloNetworks.minorreticulationgamma(net)
         R"""
@@ -210,7 +210,7 @@ function rexport(net::HybridNetwork; maintree::Bool=false, useedgelength::Bool=t
                 """
             end
         end
-    elseif net.numHybrids == 0
+    elseif net.numhybrids == 0
         R"""
         class(phy) = "phylo"
         """
diff --git a/src/substitutionmodels.jl b/src/substitutionmodels.jl
deleted file mode 100644
index d675d72..0000000
--- a/src/substitutionmodels.jl
+++ /dev/null
@@ -1,64 +0,0 @@
-plot(mod::PhyloNetworks.TraitSubstitutionModel) = error("plot not defined for $(typeof(mod)).")
-
-"""
-    plot(model::TwoBinaryTraitSubstitutionModel)
-
-Display substitution rates for a trait evolution model for two
-possibly dependent binary traits; using `R` via `RCall`.
-Adapted from fitPagel functions found in the `R` package `phytools`.
-
-# Examples
-
-```
-julia> using PhyloNetworks
-
-julia> m = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
-                ["carnivory", "noncarnivory", "wet", "dry"])
-Substitution model for 2 binary traits, with rate matrix:
-                     carnivory-wet    carnivory-dry noncarnivory-wet noncarnivory-dry
-    carnivory-wet                *           1.0000           2.0000           0.0000
-    carnivory-dry           3.1000                *           0.0000           1.1000
- noncarnivory-wet           1.2000           0.0000                *           2.0000
- noncarnivory-dry           0.0000           2.2000           1.1000                *
-
-julia> plot(m);
-```
-"""
-function plot(object::PhyloNetworks.TwoBinaryTraitSubstitutionModel)
-    R"""
-    signif<-3
-    plot.new()
-    par(mar=c(1.1,2.1,3.1,2.1))
-    plot.window(xlim=c(0,2),ylim=c(0,1),asp=1)
-    """
-    R"""
-    mtext("Two Binary Trait Substitution Model",side=3,adj=0,line=1.2,cex=1.2)
-    arrows(x0=0.15,y0=0.15,y1=0.85,lwd=2,length=0.1)
-    arrows(x0=0.2,y0=0.85,y1=0.15,lwd=2,length=0.1)
-    arrows(x0=1.6,y0=0.05,x1=0.4,lwd=2,length=0.1)
-    arrows(x0=0.4,y0=0.1,x1=1.6,lwd=2,length=0.1)
-    arrows(x0=1.8,y0=0.15,y1=0.85,lwd=2,length=0.1)
-    arrows(x0=1.85,y0=0.85,y1=0.15,lwd=2,length=0.1)
-    arrows(x0=1.6,y0=0.9,x1=0.4,lwd=2,length=0.1)
-    arrows(x0=0.4,y0=0.95,x1=1.6,lwd=2,length=0.1)
-    text(x=0.175,y=0.95,$(object.label[1]))
-    text(x=1.825,y=0.95,$(object.label[2]))
-    text(x=1.825,y=0.05,$(object.label[4]))
-    text(x=0.175,y=0.05,$(object.label[3]))
-    """
-    R"""
-    text(x=1,y=1,round($(object.rate[5]),signif),cex=0.8)
-    """
-    R"""
-    text(x=1,y=0.85,round($(object.rate[6]),signif),cex=0.8)
-    text(x=1.9,y=0.5,round($(object.rate[3]),signif),cex=0.8,srt=90)
-    text(x=1.75,y=0.5,round($(object.rate[4]),signif),cex=0.8,srt=90)
-    """
-    R"""
-    text(x=1,y=0,round($(object.rate[8]),signif),cex=0.8)
-    text(x=1,y=0.15,round($(object.rate[7]),signif),cex=0.8)
-    text(x=0.1,y=0.5,round($(object.rate[2]),signif),cex=0.8,srt=90)
-    text(x=0.25,y=0.5,round($(object.rate[1]),signif),cex=0.8,srt=90)
-    """
-    return nothing
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index 422cff4..c81cfbe 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,6 +8,5 @@ using PhyloNetworks # to read topologies
 @testset "PhyloPlots Tests" begin
   include("test_phylonetworkPlots.jl")
   include("test_plotRCall.jl")  # uses DataFrames to test annotations
-  include("test_substitutionmodels.jl")
   include("test_rexport.jl")    # uses RCall, but NOT library(ape)
 end
diff --git a/test/test_phylonetworkPlots.jl b/test/test_phylonetworkPlots.jl
index 4cc9a37..36d298c 100644
--- a/test/test_phylonetworkPlots.jl
+++ b/test/test_phylonetworkPlots.jl
@@ -1,6 +1,6 @@
 @testset "Test setup for plotting PhyloNetworks objects" begin
 
-  net = readTopology("(A:2.5,((B:1,#H1:0.5::0.1):1,(C:1,(D:0.5)#H1:0.5::0.9):1):0.5);")
+  net = readnewick("(A:2.5,((B:1,#H1:0.5::0.1):1,(C:1,(D:0.5)#H1:0.5::0.9):1):0.5);")
   # in this order:
   # edge_xB, edge_xE, edge_yB, edge_yE,
   # node_x, node_y, node_yB, node_yE,
@@ -67,7 +67,7 @@
   # example with level-2 network, non-tree child:
   # one hybrid node ends up as a leaf in the major tree.
   # no major child edge to follow to set coordinates
-  net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);")
+  net = readnewick("((((B)#H1:::0.2)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);")
   @test_logs plot(net, shownodenumber=true, showgamma=true);
   @test PhyloPlots.edgenode_coordinates(net, false, false) == (
     [5.0, 4.0, 1.0, 3.0, 3.0, 3.0, 2.0, 4.0, 4.0, 2.0, 1.0],
@@ -80,7 +80,7 @@
     [5.0, 5.0, 4.0, 2.0, 3.0, 4.0, 6.0, 6.0, 5.5, 4.25],
     [5.0, 4.0], [5.0, 4.0], [4.0, 1.0], [5.0, 4.0],
     1.0, 6.0, 1.0, 6.0)
-  net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2)S1,(#H1,A)S2)S3)S4);")
+  net = readnewick("((((B)#H1:::0.2)#H2,((D,C,#H2)S1,(#H1,A)S2)S3)S4);")
   @test_logs plot(net, shownodenumber=true, showgamma=true);
   @test PhyloPlots.edgenode_coordinates(net, false, false) == (
     [5.0, 4.0, 1.0, 3.0, 3.0, 3.0, 2.0, 4.0, 4.0, 2.0, 1.0],
diff --git a/test/test_plotRCall.jl b/test/test_plotRCall.jl
index 92fc293..c4ec63e 100644
--- a/test/test_plotRCall.jl
+++ b/test/test_plotRCall.jl
@@ -2,12 +2,10 @@
   # testing for absence of errors, not for correctness
 
   # network rooted at a leaf: test for no error in warning message
-  net = readTopology("(A:1,B:1);"); net.root = 2
+  net = readnewick("(A:1,B:1);"); net.rooti = 2
   @test_logs (:warn, r"rootonedge!\(network_name, 2\)") plot(net)
-  #net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);")
-  net = readTopology("(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);");
-  # test deprecated old function
-  @test_logs (:warn, r"is deprecated") plot(net, :R, showNodeNumber=true, showGamma=true)
+  #net = readnewick("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);")
+  net = readnewick("(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);");
 
   @test_logs plot(net);
   @test_logs plot(net, edgewidth=Dict(1=>4, 28=>28))
@@ -27,10 +25,13 @@
 
 
   # coverage for the nomajorchild
-  net2 = readTopology("((((B)#H1)#H2,((D,C,#H2)S1,(#H1:::.8,A)S2)S3)S4);")
-  @test_logs plot(net2, style=:majortree);
+  net2 = readnewick("((((B)#H1)#H2,((D,C,#H2)S1,(#H1:::.8,A)S2)S3)S4);")
+  res = (@test_logs plot(net2, style=:majortree))
+  @test keys(res) == (:xmin, :xmax, :ymin, :ymax, :node_x, :node_y,
+    :node_y_lo, :node_y_hi, :edge_x_lo, :edge_x_hi, :edge_y_lo, :edge_y_hi,
+    :node_data, :edge_data)
 
   # plot based on RCall and ape:
-  tre = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
+  tre = readnewick("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
   # fixit: plot(tre, :ape)
 end
diff --git a/test/test_rexport.jl b/test/test_rexport.jl
index c4cf5cf..a7dcb2b 100644
--- a/test/test_rexport.jl
+++ b/test/test_rexport.jl
@@ -9,8 +9,8 @@ useape = false; # these extra tests would require the installation of ape by the
 
 # on a tree with some edge lengths missing
 s = "(A,(B:1.0,(C:1.0,D:1.0):1.0):1.0);";
-tree1 = PhyloPlots.rexport(readTopology(s)); # for testing rexport
-net2 = readTopology(s); # for testing sexp on clean network
+tree1 = PhyloPlots.rexport(readnewick(s)); # for testing rexport
+net2 = readnewick(s); # for testing sexp on clean network
 if useape
 ## check for correct unrooted topology and for identical lengths
 ## of (unrooted) internal edges, but requires t
@@ -34,8 +34,8 @@ class(tree2r) = "phylo"
 
 # network, h=1, missing gamma values
 s = "(((A:4.0,(B:1.0)#H1:1.1):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);";
-phy1 = PhyloPlots.rexport(readTopology(s));
-net2 = readTopology(s);
+phy1 = PhyloPlots.rexport(readnewick(s));
+net2 = readnewick(s);
 if useape # needs ape version > 4.1 with read.evonet (not in 4.1)
     R"phyr = read.evonet(text = $s)"
     @test convert(Bool, R"dist.topo($phy1, phyr, method='score') == 0")
@@ -58,9 +58,9 @@ class(phy2r) = c("evonet","phylo")
 # network, h=1, maintree=true; minor hybrid edge length missing
 s = "(((A:4.0,(B)#H1:1.1::0.9):0.5,(C:0.6,#H1):1.0):3.0,D:5.0);";
 stree = "(((A:4.0,B):0.5,C:1.6):3.0,D:5.0);"; # main tree
-phy1 = PhyloPlots.rexport(readTopology(s); maintree=true);
-net2 = readTopology(s);
-tree2 = majorTree(readTopology(s));
+phy1 = PhyloPlots.rexport(readnewick(s); maintree=true);
+net2 = readnewick(s);
+tree2 = majortree(readnewick(s));
 if useape
     R"tree2r = read.tree(text = $stree)";
     @test convert(Bool, R"dist.topo($phy1, tree2r, method='score') == 0")
@@ -86,7 +86,7 @@ class(phy2r) = c("evonet","phylo")
 
 # on a network, h=1, with useedgelength=false
 s = "(((A:1.0,(B:.5)#H1:.2::0.9):.8,(C:1.5,#H1:.01::0.1):2):.6,D:5.0);";
-phy1 = PhyloPlots.rexport(readTopology(s), useedgelength=false);
+phy1 = PhyloPlots.rexport(readnewick(s), useedgelength=false);
 if useape # needs ape version > 4.1 with read.evonet (not in 4.1), also dist.topo for networks
     R"phyr = read.evonet(text = $s)"
     @test convert(Bool, R"dist.topo($phy1, phyr) == 0")
@@ -104,8 +104,8 @@ class(phy1r) = c("evonet","phylo")
 # on a network with h=2 hybridizations
 s = "(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,
       (((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);";
-phy1 = PhyloPlots.rexport(readTopology(s));
-net2 = readTopology(s);
+phy1 = PhyloPlots.rexport(readnewick(s));
+net2 = readnewick(s);
 R"""
 phy1r = list(Nnode=14L, edge=matrix(c(13,13,13,14,14,15,15,16,16,17,17,18,18,19,20,20,21,21,22,23,24,25,25,26,26,
              15,14,12,10,11,18,16,17,9,24,8,1,19,20,21,23,2,22,3,4,25,26,7,5,6), 25,2),
diff --git a/test/test_substitutionmodels.jl b/test/test_substitutionmodels.jl
deleted file mode 100644
index e37231f..0000000
--- a/test/test_substitutionmodels.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-@testset "test of substitution model plot, 2 binary traits" begin
-
-m3 = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
-     ["carnivory", "noncarnivory", "wet", "dry"]);
-@test_logs plot(m3)
-
-end