I know that this is a really old question, but still it comes pretty early in my google search, I though I'd give you my solution:

library(ape)
mat <- matrix(
    c(
      1,1,1,1,1,1,1,
      0,0,1,0,1,0,1,
      0,0,1,0,0,0,1,
      1,1,1,0,1,1,1,
      0,0,0,0,0,0,0

    ), byrow = TRUE, nrow = 5, dimnames = list(paste0("r", 1:5), paste0("c", 1:7))
)

mat.nj <- nj(dist.gene(mat)) # neighbour joining tree construction
plot(mat.nj, "phylo") # we plot it
nodelabels() # we id node label

# internal node state reconstruction
MPR(mat[,1], mat.nj, outgroup = "r1") # this give us the value of all internal nodes for c1 assuming root of the tree is r1

# lapply magic to call MPR on all character/gene/whatever
internalNodes <- lapply(seq_len(ncol(mat)), function(x) MPR(mat[, x], mat.nj, outgroup = "r1")[, 1, drop = FALSE])
internalNodes <- do.call(cbind, internalNodes)

The Phylip package has something. but if you find how to do it with ape, I am interested to know.