How to best summarize multivariate and unimodal but asymmetric posterior pdf generated with an MCMC?

To do parameter inference in an inverse problem, I use a Metroplois Hastings MCMC to sample the posterior pdf of the model parameters $$\theta$$ of the forward model $$f(\theta)$$. Where the forward model $$f(\cdot)$$ maps a vector $$\theta$$ in the parameter space $$X$$ to a vector $$\nu$$ in the space of observable features $$Y$$.

When sampling the pdf $$p(\theta|\nu_\text{obs})$$ using the MCMC, where $$\nu_\text{obs}$$ is a noisy observation in $$Y$$, I get a multivariate and unimodal pdf. Though the marginalized pdf of $$p(\theta|\nu_\text{obs})$$ for some of the elements of $$\theta$$ are highly asymmetric.

How do I best summarize such a posterior? For symmetric posteriors I would usually use the median an the 16th and 84th percentiles. For univariate asymmetric pdf I usually use the highest density interval (HDI) and the mode.

I know that there are packages to calculate the HDI of bivariate distributions, but I am not aware of any for more dimensions.