The construction is functorial with respect to algebraic morphisms of Lie algebroids (as opposed to geometric ones): see for instance my paper with Van de Bergh https://arxiv.org/pdf/0708.2725.pdf (it encompasses the globalization techniques that are used in the paper you cite).
As Stefan pointed out, the choice of a connection is rather important for the Dolgushev-Fedosov resolutions. In the paper with Van den Bergh we use much bigger resolutions that are independant of such a choice (they are somehow universal). A torsion free connection allows one to linearize the jet bundle. The resolution we construct in the paper with Van den Bergh involves an algebra that is universal among the ones that can linearize the jet bundle (and as such it somehow carries a universal connection). This in particular shows that the homotopy class of the $L_\infty$-morphism that you get from my paper does not depend on the choice of the connection involved in the construction.
Note that the algebraic morphisms from the paper with Van den Bergh are called comorphisms by several people (see e.g. https://arxiv.org/pdf/1210.4443.pdf).
Let me also observe that if you fix the base and only consider morphisms that are the identity on the base, then the categoy of Lie algebroids with comorphisms/algebraic morphisms is just the opposite category to the category of Lie algebroids with morphisms/geometric morphisms.
So, combining all these comments, I would say that the answer to your question is yes.
EDIT JULY 25 2017: the functoriality only works for isomorphisms of Lie algebroids (see comments below).
- the functoriality only works for isomorphisms of Lie algebroids (see comments below), and more generally for étale morphisms (like open embeddings in the tangent Lie algebroid case).
- there are probably other situations where it works. As suggested in the question, one could expect it to work for inclusions of Lie algebroids.