I believe that the answer to your question about whether or not connections are known in higher dimensions is yes. Before I get into the details however, it might be useful to recall what happens for quaternion algebras.

Let $k$ be a totally real field and $B$ be a quaternion division algebra over $k$ which is split at a unique real place of $k$. This split place, call it $\mathfrak p$, furnishes us with an isomorphism $B_{\mathfrak p} \cong M_2(\mathbf{R})$. Let $\mathcal O$ be a maximal order of $B$ and $\mathcal O^1$ be the multiplicative subgroup consisting of elements with reduced norm $1$. Let $\Gamma_\mathcal O$ denote the image in $PSL_2(\mathbf{R})$ of $\mathcal O^1$. Recalling that $PSL_2(\mathbf{R})$ is the group of orientation-preserving isometries of the hyperbolic plane, the work of Borel and Harish-Chandra imply that the quotient $\mathbf{H}^2/\Gamma_\mathcal O$ is a compact hyperbolic $2$-orbifold of finite area.

Note that everything in the last paragraph works equally well for other number fields. If $k$ were a number field with a unique complex place and $B$ was ramified at all real places of $k$, then we would end up getting lattices in $PSL_2(\mathbf{C})$, hence hyperbolic $3$-orbifolds. More generally one gets lattices in products of $PSL_2(\mathbf{R})$ and $PSL_2(\mathbf{C})$. These manifolds (or orbifolds) will only be hyperbolic in the two special cases I mentioned.

Everything that I have written thus far generalizes to arbitrary dimensions in the obvious way to yield lattices in groups like $PSL_d(\mathbf{R})$, $PSL_d(\mathbf{C})$ and their products. The manifolds (or orbifolds) one gets will therefore be quotients of the symmetric spaces of these groups.

You also asked about the spectral theory of these higher dimensional manifolds and the embedding theory of orders in the relevant division algebra. In this context the selectivity theory you would want to use is due to Luis Arenas-Carmona (http://arxiv.org/abs/1403.5826). Geometrically however, things are a bit trickier.

I think that the first person to seriously do anything with the spectral theory of orbifolds arising from these sort of constructions was Marie-France Vignéras (Variétés riemanniennes isospectrales et non isométriques, Ann. of Math. (2) 112 (1980), no. 1, 21–32. ) In this paper Vignéras constructed arbitrarily large families of isospectral but not isometric manifolds with universal cover the symmetric space of products of $PSL_2(\mathbf{R})$ and $PSL_2(\mathbf{C})$. Note that this stuff is also written up in Chapter 12 of Maclachlan and Reid's book "The arithmetic of hyperbolic 3-manifolds". The point is that what is actually going on behind the scenes is that a trace formula is being used. In this case, it follows from the Selberg Trace Formula that two hyperbolic surfaces are isospectral (with respect to the Laplace operator) if and only if they have the same geodesic length spectra (e.g. sets of lengths of closed geodesics). It is the latter that is intimately connected with orders in quaternion algebras, for geodesic lengths on the orbifold obtained from $\Gamma_\mathcal O $ correspond to conjugacy classes of elements in $\Gamma_\mathcal O$ with a given minimal polynomial, which in turn correspond to conjugacy classes of embeddings of certain quadratic $\mathcal O_k$-orders into $\mathcal O$. Thus for $\Gamma_\mathcal O$ and $\Gamma_{\mathcal O'}$ to yield isospectral orbifold quotients, one would need to know that $\mathcal O$ and $\mathcal O'$ admit embeddings of precisely the same set of quadratic orders. And this of course is what the notion of selectivity in quaternion algebras was designed to address.

If you wanted to use this sort of method to produce isospectral orbifolds from higher dimensional division algebras and other variants thereof (for instance, central simple algebras equipped with certain types of involutions), you would need a trace formula to tell you that it sufficed to produce orbifolds with the same geodesic length spectra. Unfortunately these trace formulas are only known in a small number of cases and there are definitely cases in which the Laplace spectrum of a (locally-symmetric) manifold determines the lengths which occur in the geodesic length spectrum but not their multiplicities.

So you would probably want to study geodesics on your manifolds. This is still doable in this more general setting, but more difficult. Also keep in mind the paper "Division algebras and non-commensurable isospectral manifolds" (http://arxiv.org/abs/math/0501064) by Lubotzky, Samuels and Vishne. In this paper they construct isospectral manifolds from higher dimensional division algebras. In fact they manage to do so without studying orders and embeddings at all. The idea is based upon the fact that a division algebra and its opposite algebra will have the exactly the same maximal subfields but will not be isomorphic when the dimension is greater than $4$. The authors consider lattices in a division algebra and its opposite algebra, and then use the global Jacquet-Langlands correspondence to deduce that their manifolds have equal Laplace spectra. Note that when you use the Vigneras style construction of isospectral manifolds you always get manifolds that are commensurable, hence "less interesting" than the Lubotzky, Samuels and Vishne examples.

Of course there are still plenty of very interesting things that one can say about the geometry of the manifolds obtained from orders in higher dimensional division algebras. No references come to mind immediately, though for some inspiration you might try looking at:

- Chapter 10 of Maclachlan and Reid
- The papers of Prasad and Rapinchuk. You might start with: http://arxiv.org/abs/0809.2401 and the references therein. Note that Prasad and Rapinchuk are experts on algebraic groups and so work in a lot more generality than you would want to. If you are not an expert in this style of thinking, it might be a good exercise to try to work out exactly what happens when the groups being considered are obtained from division algebras defined over a number field.