Although HCubature.jl does indeed support arbitrary precision (its multidimensional quadrature rule is carefully written to use the precision of the endpoints), it corresponds to a low-order quadrature rule that will be very slow (require a lot of subdivisions) to reach high accuracies.
If you are using BigFloat or similar types to obtain high accuracy (\gg 16 digits), then you really want to use a high-order quadrature rule (computed in BigFloat precision too). The QuadGK.jl package supports this, as discussed in this tutorial example and this post, for example. You can use nested quadgk calls to do multidimensional integration, and there is a package IteratedIntegration.jl wrapping QuadGK that does this more efficiently.
However, if you want high accuracy arbitrary precision integration, you should also have an integrand that is very smooth (or for which any singularities are built into the quadrature rule). In that case, you might not want an adaptive routine at all, but simply use a tensor product of 1d quadrature rules (computed with QuadGK to arbitrary precision), e.g. see the example in this post.