This is mostly answered by decision theory. The approach is that you must formalize what you want to achieve with your point estimate, and that can be formalized as a loss function. Examples are squared error loss $$ L(\hat{\theta}, \theta) = (\hat{\theta} - \theta)^2, $$ absolute value loss, or others. Then the optimal estimator (in your Bayesian setting) is given by the function $\hat{\theta}$ that minimized posterior expected loss $$ E L(\hat{\theta}, \theta) $$ where the expectation is over the posterior distribution of $\theta$. In the case of the squared error loss, that is the posterior expected value.
You can read about this an any text about Bayesian decision theory, such as Christian Robert: "The Bayesian choice" (Springer).
EDIT
As for the question/commentary in the comment: The loss function is supposed to give your loss structure, and no paper can tell you that! Of course, if yours is an inference problem, so there is no immediate economical or other loss to you, the answer is not easy, but the book I referred to above will contain some discussion. As for the blog you referred, http://doingbayesiandataanalysis.blogspot.se/2012/04/why-to-use-highest-density-intervals.html that is an interestingt point. It advises for highest posterior density intervals (HPI). As far as I understand, that will not follow as conclusion from any decision-theoretic approach. HPI-intervals are not invariant under reparametrization, any proper Bayesian solution must be invariant. The meaning of this is that if your parameter is $\theta$, and you reparametrize the model with $g(\theta)$, where $g$ is some invertible function, then you get the same result directly doing inference for $\theta$, or by first reformulating everything in terms of $g(\theta)$, doing inference on that new scale, and then using $g^{-1}$ to translate that inference back to the original scale $\theta$. That will not work out for HPI, so is not given by any loss function.
