diff --git a/.github/workflows/R-CMD-check.yml b/.github/workflows/R-CMD-check.yml
index ec2d37fa6..8a4e128d5 100644
--- a/.github/workflows/R-CMD-check.yml
+++ b/.github/workflows/R-CMD-check.yml
@@ -84,10 +84,10 @@ jobs:
           remotes::install_cran("rcmdcheck")
         shell: Rscript {0}
         
-      - name: Install vdiffr for checks (macOS)
+      - name: Install covr, vdiffr for checks (macOS)
         if: runner.os == 'macOS'
         run: |
-          remotes::install_cran("vdiffr", type = 'source')
+          remotes::install_cran(c("covr", "vdiffr"))
           library('vdiffr')
         shell: Rscript {0}
 
diff --git a/DESCRIPTION b/DESCRIPTION
index b57db43fb..157e8b466 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: TreeSearch
 Title: Phylogenetic Tree Search Using Custom Optimality Criteria
-Version: 0.4.3.9008
+Version: 0.4.3.9009
 Authors@R: c(person("Martin R.", 'Smith',
                   email="martin.smith@durham.ac.uk", 
                   role=c("aut", "cre", "cph"), 
diff --git a/NAMESPACE b/NAMESPACE
index 0e9d09132..88c3f32e7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -19,6 +19,7 @@ export(ClusteringConcordance)
 export(ConcordantInfo)
 export(ConcordantInformation)
 export(Consistency)
+export(ConstrainedNJ)
 export(DoNothing)
 export(EasyTrees)
 export(EasyTreesy)
diff --git a/NEWS.md b/NEWS.md
index 248d33ff2..48e90df04 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,10 +1,14 @@
-# TreeSearch 0.4.3.9008 (development)
+# TreeSearch 0.4.3.9009 (development)
+
+## Changes in functionality
 
  * v......9007 introduces a pre-ratchet TBR search and adjusts default settings.
  * v0.4.3.9005 fixes some bugs that may affect tree scoring under IW & PP.
 
+## New functions
  - `EasyTrees()` 'shiny' graphical user interface for tree search
- - `PlotCharacter()` reconstructs character distributions on trees.
+ - `PlotCharacter()` reconstructs character distributions on trees
+ - `ConstrainedNJ()` constructs starting trees that respect a constraint
  - `SiteConcordance()` calculates exact site concordance
  - `ConcordantInformation()` evaluates signal:noise of dataset implied by a
    given tree.
@@ -14,7 +18,9 @@
  - `Consistency()` calculates consistency and retention 'indices'
  - `MinimumLength()` and `MaximumLength()` calculate range of possible lengths
    of characters in a dataset on any tree
- - `TreeLength()` supports lists of trees.
+
+## Improvements
+ - `TreeLength()` supports lists of trees
  - Set handling of 'gap' token (-) when creating Morphy object with `gap = `
  - Label nodes with split frequencies using `JackLabels(plot = FALSE)`
  - Support for topological constraints during tree search
diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R
index 76fdd9fec..51809a21f 100644
--- a/R/MaximizeParsimony.R
+++ b/R/MaximizeParsimony.R
@@ -30,7 +30,8 @@
 #' @param tree (optional) A bifurcating tree of class \code{\link{phylo}},
 #' containing only the tips listed in `dataset`, from which the search
 #' should begin.
-#' If unspecified, a neighbour-joining tree will be generated from `dataset`.
+#' If unspecified, a neighbour-joining tree will be generated from `dataset`,
+#' respecting any supplied `constraint`.
 #' Edge lengths are not supported and will be deleted.
 #' @param ratchIter Numeric specifying number of iterations of the 
 #' parsimony ratchet (Nixon 1999) to conduct.
@@ -61,9 +62,8 @@
 #' Setting to larger values will include trees suboptimal by up to `tolerance`
 #' in search results, which may improve the accuracy of the consensus tree
 #' (at the expense of resolution) (Smith 2019).
-#' @param constraint Either `NULL` or an object of class `phyDat`. Trees that
-#' are not perfectly compatible with each character in `constraint` will not
-#' be considered during search.
+#' @param constraint An object of class `phyDat`; returned trees will be
+#' perfectly compatible with each character in `constraint`.
 #' See [vignette](https://ms609.github.io/TreeSearch/articles/inapplicable.html)
 #' for further examples.
 #' @param verbosity Integer specifying level of messaging; higher values give
@@ -157,7 +157,7 @@
 #' \insertRef{Smith2019}{TreeSearch}
 #' @encoding UTF-8
 #' @export
-MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
+MaximizeParsimony <- function (dataset, tree,
                                ratchIter = 6L,
                                tbrIter = 2,
                                startIter = 2L, finalIter = 1L,
@@ -166,7 +166,7 @@ MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
                                quickHits = 1 / 3,
                                concavity = Inf,
                                tolerance = sqrt(.Machine$double.eps),
-                               constraint = NULL,
+                               constraint,
                                verbosity = 2L, session = NULL) {
   # Define functions
   .Message <- if (is.null(session)) function (level, ...) {
@@ -463,8 +463,23 @@ MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
     class = 'multiPhylo')
   }
   
+  
+  # Define constants
+  epsilon <- tolerance #sqrt(.Machine$double.eps)
+  profile <- .UseProfile(concavity)
+  iw <- is.finite(concavity)
+  constrained <- !missing(constraint)
+  startTime <- Sys.time()
+  stopTime <- startTime + as.difftime(maxTime, units = 'mins')
+  
   # Initialize tree
-  if (inherits(tree, 'multiPhylo')) {
+  if (missing(tree)) {
+    if (constrained) {
+      tree <- ConstrainedNJ(dataset, constraint)
+    } else {
+      tree <- NJTree(dataset)
+    }
+  } else if (inherits(tree, 'multiPhylo')) {
     .Message(1L, "Starting search from `tree[[1]]`.")
     tree <- tree[[1]]
   }
@@ -487,13 +502,6 @@ MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
     outgroup <- NA
   }
   
-  # Define constants
-  epsilon <- tolerance #sqrt(.Machine$double.eps)
-  profile <- .UseProfile(concavity)
-  iw <- is.finite(concavity)
-  constrained <- !is.null(constraint)
-  startTime <- Sys.time()
-  stopTime <- startTime + as.difftime(maxTime, units = 'mins')
   
   # Initialize constraints
   if (constrained) {
@@ -859,7 +867,7 @@ Resample <- function (dataset, tree = NJTree(dataset), method = 'jack',
                       ratchIter = 1L, tbrIter = 8L, finalIter = 3L,
                       maxHits = 12L, concavity = Inf,
                       tolerance = sqrt(.Machine$double.eps),
-                      constraint = NULL,
+                      constraint,
                       verbosity = 2L, session = NULL,
                       ...) {
   if (!inherits(dataset, 'phyDat')) {
@@ -897,6 +905,33 @@ Resample <- function (dataset, tree = NJTree(dataset), method = 'jack',
                     verbosity = verbosity, session = session, ...) 
 }
 
+#' Constrained neighbour-joining tree
+#' 
+#' Constructs a neighbour-joining tree such that the tree is consistent with a
+#' constraint.
+#' 
+#' @param weight Numeric specifying degree to upweight characters in
+#' `constraint`.
+#' 
+#' @return `ConstrainedNJ()` returns a tree of class `phylo`.
+#' @importFrom TreeTools NJTree
+#' @inheritParams MaximizeParsimony
+#' @examples
+#' dataset <- TreeTools::MatrixToPhyDat(matrix(
+#'   c(0, 1, 1, 1, 0, 1,
+#'     0, 1, 1, 0, 0, 1), ncol = 2,
+#'   dimnames = list(letters[1:6], NULL)))
+#' constraint <- TreeTools::MatrixToPhyDat(
+#'   c(a = 0, b = 0, c = 0, d = 0, e = 1, f = 1))
+#' plot(ConstrainedNJ(dataset, constraint))
+#' @template MRS
+#' @export
+ConstrainedNJ <- function (dataset, constraint, weight = 12345) {
+  conData <- c(constraint, dataset)
+  attr(conData, 'weight')[seq_len(attr(constraint, 'nr'))] <- weight
+  NJTree(conData)
+}
+
 #' Launch tree search graphical user interface
 #' 
 #' @rdname MaximizeParsimony
diff --git a/man/ConstrainedNJ.Rd b/man/ConstrainedNJ.Rd
new file mode 100644
index 000000000..1393fd8ad
--- /dev/null
+++ b/man/ConstrainedNJ.Rd
@@ -0,0 +1,40 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/MaximizeParsimony.R
+\name{ConstrainedNJ}
+\alias{ConstrainedNJ}
+\title{Constrained neighbour-joining tree}
+\usage{
+ConstrainedNJ(dataset, constraint, weight = 12345)
+}
+\arguments{
+\item{dataset}{A phylogenetic data matrix of class \code{\link[phangorn]{phyDat}},
+whose names correspond to the labels of any accompanying tree.}
+
+\item{constraint}{An object of class \code{phyDat}; returned trees will be
+perfectly compatible with each character in \code{constraint}.
+See \href{https://ms609.github.io/TreeSearch/articles/inapplicable.html}{vignette}
+for further examples.}
+
+\item{weight}{Numeric specifying degree to upweight characters in
+\code{constraint}.}
+}
+\value{
+\code{ConstrainedNJ()} returns a tree of class \code{phylo}.
+}
+\description{
+Constructs a neighbour-joining tree such that the tree is consistent with a
+constraint.
+}
+\examples{
+dataset <- TreeTools::MatrixToPhyDat(matrix(
+  c(0, 1, 1, 1, 0, 1,
+    0, 1, 1, 0, 0, 1), ncol = 2,
+  dimnames = list(letters[1:6], NULL)))
+constraint <- TreeTools::MatrixToPhyDat(
+  c(a = 0, b = 0, c = 0, d = 0, e = 1, f = 1))
+plot(ConstrainedNJ(dataset, constraint))
+}
+\author{
+\href{https://smithlabdurham.github.io/}{Martin R. Smith}
+(\href{mailto:martin.smith@durham.ac.uk}{martin.smith@durham.ac.uk})
+}
diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd
index bc52cafb0..2e0dab159 100644
--- a/man/MaximizeParsimony.Rd
+++ b/man/MaximizeParsimony.Rd
@@ -10,7 +10,7 @@
 \usage{
 MaximizeParsimony(
   dataset,
-  tree = NJTree(dataset),
+  tree,
   ratchIter = 6L,
   tbrIter = 2,
   startIter = 2L,
@@ -20,7 +20,7 @@ MaximizeParsimony(
   quickHits = 1/3,
   concavity = Inf,
   tolerance = sqrt(.Machine$double.eps),
-  constraint = NULL,
+  constraint,
   verbosity = 2L,
   session = NULL
 )
@@ -36,7 +36,7 @@ Resample(
   maxHits = 12L,
   concavity = Inf,
   tolerance = sqrt(.Machine$double.eps),
-  constraint = NULL,
+  constraint,
   verbosity = 2L,
   session = NULL,
   ...
@@ -53,7 +53,8 @@ whose names correspond to the labels of any accompanying tree.}
 \item{tree}{(optional) A bifurcating tree of class \code{\link{phylo}},
 containing only the tips listed in \code{dataset}, from which the search
 should begin.
-If unspecified, a neighbour-joining tree will be generated from \code{dataset}.
+If unspecified, a neighbour-joining tree will be generated from \code{dataset},
+respecting any supplied \code{constraint}.
 Edge lengths are not supported and will be deleted.}
 
 \item{ratchIter}{Numeric specifying number of iterations of the
@@ -94,9 +95,8 @@ Setting to larger values will include trees suboptimal by up to \code{tolerance}
 in search results, which may improve the accuracy of the consensus tree
 (at the expense of resolution) (Smith 2019).}
 
-\item{constraint}{Either \code{NULL} or an object of class \code{phyDat}. Trees that
-are not perfectly compatible with each character in \code{constraint} will not
-be considered during search.
+\item{constraint}{An object of class \code{phyDat}; returned trees will be
+perfectly compatible with each character in \code{constraint}.
 See \href{https://ms609.github.io/TreeSearch/articles/inapplicable.html}{vignette}
 for further examples.}
 
diff --git a/tests/testthat/test-MaximizeParsimony.R b/tests/testthat/test-MaximizeParsimony.R
index 7b01abcf6..0e0b21140 100644
--- a/tests/testthat/test-MaximizeParsimony.R
+++ b/tests/testthat/test-MaximizeParsimony.R
@@ -35,6 +35,17 @@ test_that("Constraints work", {
   
 })
 
+test_that("Constrained NJ trees work", {
+  dataset <- MatrixToPhyDat(matrix(
+   c(0, 1, 1, 1, 0, 1,
+     0, 1, 1, 0, 0, 1), ncol = 2,
+   dimnames = list(letters[1:6], NULL)))
+  constraint <- MatrixToPhyDat(c(a = 0, b = 0, c = 0, d = 0, e = 1, f = 1))
+  expect_equal(ape::read.tree(text = "(a, (d, ((c, b), (e, f))));"),
+               ConstrainedNJ(dataset, constraint))
+  expect_equal(NJTree(dataset), ConstrainedNJ(dataset, dataset))
+})
+
 test_that("Inconsistent constraints fail", {
   constraint <- MatrixToPhyDat(matrix(
     c(0, 1, 1, 1, 0, 0,