diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml
index 8c68f5464..e02d360ca 100644
--- a/.github/workflows/coverage.yml
+++ b/.github/workflows/coverage.yml
@@ -54,6 +54,9 @@ jobs:
           needs: |
             coverage
 
+      - name: Set up pandoc
+        uses: r-lib/actions/setup-pandoc@v2
+
       - name: Code coverage
         run: |
           covr::package_coverage(quiet = FALSE) # https://github.com/r-lib/covr/issues/252
diff --git a/DESCRIPTION b/DESCRIPTION
index f688aa6b6..dd8cc1a7d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: TreeSearch
 Title: Phylogenetic Analysis with Discrete Character Data
-Version: 1.5.1.9002
+Version: 1.5.1.9003
 Authors@R: c(
     person(
       "Martin R.", 'Smith',
diff --git a/NEWS.md b/NEWS.md
index 471dfabfc..333889154 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# TreeSearch 1.5.1.9003 (2025-02)
+
+- Improve support for constraints in `AdditionTree()`
+  [#173](https://github.com/ms609/TreeSearch/issues/173)
+
 # TreeSearch 1.5.1.9002 (2025-01)
 
 - Buttons to download consensus trees in app
diff --git a/R/AdditionTree.R b/R/AdditionTree.R
index 29139ec1d..d8cb340d4 100644
--- a/R/AdditionTree.R
+++ b/R/AdditionTree.R
@@ -31,7 +31,7 @@ AdditionTree <- function (dataset, concavity = Inf, constraint, sequence) {
   # Initialize missing parameters
   taxa <- names(dataset)
   if (missing(sequence)) {
-    sequence <- taxa[1]
+    sequence <- taxa[[1]]
   } else if (is.numeric(sequence)) {
     sequence <- taxa[sequence]
   }
@@ -78,7 +78,7 @@ AdditionTree <- function (dataset, concavity = Inf, constraint, sequence) {
         constraint <- MatrixToPhyDat(t(as.matrix(constraint)))
       }
       thisConstr <- constraint[theseTaxa]
-      if (length(thisConstr[[1]]) && min(table(unlist(thisConstr))) > 1) {
+      if (.ConstraintConstrains(thisConstr)) {
         # Constraint constrains theseTaxa
         
         morphyConstr <- PhyDat2Morphy(thisConstr)
@@ -114,6 +114,31 @@ AdditionTree <- function (dataset, concavity = Inf, constraint, sequence) {
   tree
 }
 
-.Recompress <- function (dataset) {
+
+.ConstraintConstrains <- function(constraint) {
+  if (length(constraint[[1]]) < 1) {
+    FALSE
+  } else {
+    contrast <- attr(constraint, "contrast")
+    if (dim(contrast)[[2]] < 2) {
+      FALSE
+    } else {
+      cont <- `mode<-`(contrast, "logical")
+      nLevel <- dim(contrast)[[1]]
+      # Could be > 2× more efficient using lower.tri
+      exclude <- vapply(seq_len(nLevel), function(i) {
+        colSums(apply(cont, 1, `&`, cont[i, ])) == 0
+      }, logical(nLevel))
+      
+      # TODO Validate; passes existing tests, but these do not include all 
+      # edge cases, e.g. 02 03 1 1
+      splits <- exclude * tabulate(unlist(constraint), nLevel)
+      any(splits[lower.tri(splits)] > 1 & t(splits)[lower.tri(splits)] > 1)
+    }
+  }
+}
+
+
+.Recompress <- function(dataset) {
   MatrixToPhyDat(PhyDatToMatrix(dataset))
 }
diff --git a/tests/testthat/test-AdditionTree.R b/tests/testthat/test-AdditionTree.R
index 3c943504a..395e472cf 100644
--- a/tests/testthat/test-AdditionTree.R
+++ b/tests/testthat/test-AdditionTree.R
@@ -20,6 +20,50 @@ test_that("Addition tree is more parsimonious", {
   expect_lt(Score(pr, "pr"), Score(nj10, "pr"))
 })
 
+test_that(".ConstraintConstrains() succeeds", {
+  expect_false(.ConstraintConstrains(NULL))
+  
+  # Single level
+  expect_false(.ConstraintConstrains(
+    structure(list(A = 1L, B = 2L, C = 2L, D = 2L), weight = 1L, nr = 1L,
+              nc = 1L, index = 1L, levels = 0, allLevels = c("0", "?"),
+              type = "USER", contrast = 
+                structure(c(1, 1), dim = 2:1, dimnames = list(NULL, 0)),
+              class = "phyDat")
+  ))
+  
+  expect_false(.ConstraintConstrains(
+    structure(list(A = 1L, B = 2L, C = 1L, D = 1L, E = 3L), weight = 1L, nr = 1L,
+              nc = 2L, index = 1L, levels = 0:1,
+              allLevels = c("0", "1", "?"), type = "USER",
+              contrast = structure(c(1, 0, 1, 0, 1, 1), dim = 3:2,
+                                   dimnames = list(NULL, 0:1)),
+              class = "phyDat")
+  ))
+  expect_true(.ConstraintConstrains(structure(
+    list(A = 1L, B = 2L, C = 1L, D = 1L, E = 3L, F = 2L), weight = 1L, nr = 1L,
+    nc = 2L, index = 1L, levels = 0:1, allLevels = c("0", "1", "?"),
+    type = "USER", contrast = structure(c(1, 0, 1, 0, 1, 1), dim = 3:2,
+                                        dimnames = list(NULL, 0:1)),
+    class = "phyDat")
+  ))
+  expect_false(.ConstraintConstrains(structure(
+    list(A = 1L, B = 2L, C = 1L, D = 1L, E = 3L, F = 2L), weight = 1L, nr = 1L,
+    nc = 2L, index = 1L,
+    levels = 0:2, allLevels = c("0", "1", "2", "?"), type = "USER",
+    contrast = structure(c(1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1),
+                         dim = c(4, 3), dimnames = list(NULL, 0:2)),
+    class = "phyDat")
+  ))
+  expect_true(.ConstraintConstrains(structure(
+    list(A = 1L, B = 2L, C = 1L, D = 1L, E = 3L, F = 2L), weight = 1L, nr = 1L,
+    nc = 2L, index = 1L, levels = 0:2, allLevels = c("0", "1", "2", "?"),
+    type = "USER", contrast = structure(c(1, 0, 1, 1, 0, 1, 0, 1, 1),
+                                        dim = c(3, 3), dimnames = list(NULL, 0:2)),
+    class = "phyDat")
+  ))
+})
+
 test_that("Addition tree obeys constraints", {
   dataset <- MatrixToPhyDat(matrix(
     c(0, 1, 1, 1, 0, 1,