In this paper, we consider the following slightly subcritical Choquard equation \(\begin{aligned} {\left\{ \begin{array}{ll} -\Delta u=\left( \int _{\Omega }\frac{u^{2^*_{\alpha }-\varepsilon }(y)}{|x-y|^\alpha }dy\right) u^{2^*_{\alpha }-1-\varepsilon },\quad u>0,\ \ & \text{ in }\ \Omega ,\\ \quad \ \ u=0, \ \ & \text{ on }\ \partial \Omega , \end{array}\right. } \end{aligned}\) where \(N\ge 3\) , \(\Omega \) is a smooth bounded domain in \(\mathbb {R}^{N}\) , \(\alpha \in (0,N)\) , \(2^*_{\alpha }:=\frac{2N-\alpha }{N-2}\) is the upper critical exponent in the sense of the Hardy-Littlewood-Sobolev inequality, and \(\varepsilon >0\) is a small parameter. For the energy-minimizing solutions, we characterize their asymptotic profiles, derive their blow-up rates, and prove that their blow-up points are critical points of the Robin function. Moreover, through orthogonal decomposition, we obtain refined estimates for both the blow-up rate and concentration speed. For the least energy solutions, which form a particular family of the energy-minimizing solutions, we establish a precise energy expansion and further prove that the blow-up point globally minimizes the Robin function. The main results of this paper demonstrate that the asymptotic behavior of solutions to the above equation depends on both the spatial dimension N and the geometry of domain \(\Omega \) . In contrast to the local case studied in the literature, the nonlocal Choquard term brings several difficulties, and new estimates are required.