If you model a box as a point $(b_1,b_2, \ldots, b_d)$, and you define the dominance relationship $p \prec q$ $\iff$ $p_i < q_i$, for all $i$, then you are looking for the longest chain in this partial order. In $d$ dimensions this problem can be solved in $O(n \log^{d} n)$ times, I believe.
I would sketch a somewhat inferior solution, because it is late.
To see how to do it, consider the two dimensional case. For a point $p$, let $||p||_1 = |x(p)| + |y(p)|$. Consider the input points, $p_1, \ldots, p_n$, and observe that $p_i \prec p_j$ only if $||p_i||_1 < ||p_j||_1$. As such, sort the points by their $L_1$ norm, in increasing ordering, and let $p_1,\ldots, p_i$ be the points handled so far, where for each $p_j$, $j=1,\ldots, i$, we computed $\ell(p_j)$ -- the length of the longest path in the relation that ends at $p_j$. We add the extra dimension $\ell(p_j)$ to the point $p_j$, so that this is now its third coordinate.
Now, given a new point $p_{i+1}$, we need to compute the highest point (in the third dimension) that it dominates (in the first two dimensions) -- this is an open quadrant question, and can be solved in $O( \log^{O(1)} n)$ time using orthogonal range searching techniques. Assume you get the point $p_k$ as returned by the query -- we have that $\ell_{i+1} = \ell(p_{i+1}) = \ell(p_k)+1$. We need to insert this new point $(p_{i+1}, \ell_{i+1})$ in the the orthogonal range searching data-structure.
Thus, you can handle each point in $O(\log^{O(1)} n)$ time which readily implies that you can compute the longest chain in $O( n\log^{O(1)} n)$ time. Being more careful about the details, the running time can be improved to $O(n \log n)$ in two dimensions.
I have a paper on a related problem here:
https://arxiv.org/abs/1910.07586
I think a solution for your problem should be readily implied by some previous work.
More maybe later...
