Perhaps this is trivial and not very useful but if you run NUTS and come up with samples of both \theta and Z, then if you simply ignore the samples of Z then the resulting samples of \theta will be drawn from the marginal distribution p(\theta \vert D) and will not be specific to any single value of Z. This marginal distribution will be affected by choice of priors via the integration you listed however, so using a N(0,10) prior will give differing results from N(0,1).
Also, note that sampling Z, fixing its value in the model, and running NUTS for each of these sampled values should give you the same results as just running NUTS once. If you do the former, you’re mixing MCMC and vanilla Monte Carlo but you won’t be changing anything fundamental about the target distribution. Edit - not right.