There are - even after clarifications in the comments - several ways to interpret the question. So let me start out with the way that fits all parameters here:

Suppose we have a Hilbert space $\mathcal{H}$ (for qubits, $\mathcal{H}=\mathbb{C}^2$) and suppose we take a (pure) state in the Hilbert space $\mathcal{H}\otimes \mathcal{H}$.

Question 1: If we take qubits, are there only four entangled states? The answer is: Most definitely not. There are infinitely many entangled states in any dimension. Your are talking about Bell states. These are maximally entangled states, but even if we restrict your question by asking "Are there only four maximally entangled qubit states?", the answer is still "No, there are infinitely many maximally entangled states". We will see later that all maximally entangled states are connected by a unitary - and since there are infinitely many unitaries, there are infinitely many states.

Question 2: Can there be states such that the reduced density matrices are different? The answer is: Yes and no. I can easily give you an example where the reduced density matrices look the same, but the two matrices have to diagonally equivalent (their eigenvalues have to be the same including multiplicities). At this point, another one of the assumptions that where not clearly stated enters: The state $|\psi\rangle$ was supposed to be pure. For mixed states, the story is completely different.

Question 3: Is this phenomenon dependent on the dimension? No, it is not. However, there is one thing I may point out, namely that if it were true in some dimension, it would most definitely be true in all higher dimensions due to the fact that you can "just not use" the other dimensions.

Below, I will write down a state $|\psi_{sasha}\rangle$ that has all the features you ask for.

Now this was all "blabla". Let me add a bit mathematics to it, but since all these questions (at least for pure states) boil down to a single observation on entangled states, I wanted to first separately answer your questions.

This is the observation we need:

Lemma: Let $\mathcal{H}$ be a $d$-dimensional Hilbert space and let $\psi \in \mathcal{H}\otimes \mathcal{H}$. We fix a maximally entangled state $|\Omega\rangle:=\frac{1}{\sqrt{d}}\sum_{i=1}^d |i\rangle|i\rangle$ where $|i\rangle$ is an orthonormal basis of the Hilbert space. Then there exists a linear map $R:\mathcal{H}\to\mathcal{H}$ such that

$$ |\psi\rangle = (1\otimes R)|\Omega\rangle $$

Proof: The idea is: We can write $|\psi\rangle=\sum_{i=1}^k \sqrt{\lambda_i} |e_i\rangle |f_i\rangle$ for some orthonormal bases $|e_i\rangle$, $|f_i\rangle$, some $k\leq n$ and some $\lambda_i\geq 0$. This is true by virtue of the Schmidt decomposition. Since for any unitary, $U\otimes U|\Omega\rangle =|\Omega\rangle$, we can suppose that $|e_i\rangle=|i\rangle$.

Then write: $$(1\otimes R)|\Omega\rangle = \sum_{i=1}^n \frac{1}{\sqrt{d}} |i\rangle\otimes R|i\rangle \stackrel{!}{=} \sum_{i=1}^k \sqrt{\lambda_i} |i\rangle\otimes |f_i\rangle$$ When you compare the two sides, it is clear that setting $R|i\rangle:=\sqrt{\lambda_i\cdot d}|f_i\rangle$ for $i\leq k$ and $R|i\rangle=0$ for $i>k$ defines a valid linear map that proves our theorem. End of proof

This is all we need to answer your question, but for convenience, let's calculate the reduced density matrices:

Lemma: Setting as above, $|\psi\rangle=\sum_{i=1}^k \sqrt{\lambda_i}|e_i\rangle|f_i\rangle$. Then the reduced density matrices are:

$$ \rho_1(|\psi\rangle\langle \psi|)=\sum_{i=1}^k \lambda_i |e_i\rangle\langle e_i| \\ \rho_2(|\psi\rangle\langle \psi|)=\sum_{i=1}^k \lambda_i |f_i\rangle\langle f_i| $$

I'm not going to do the proof, because it is only really easy calculations. But this immediately answers your second and third question: Can there be pure states with different reduced density matrices in dimension two or higher? The answer is "yes", since you just have to choose $|e_i\rangle\neq |f_i\rangle$ for some $i$. Now, the two reduced density matrices are in their spectral decomposition (the $e_i$ and $f_i$ are actually the eigenvectors due to the fact that they are an orthogonal basis, which means that the $\lambda_i$ are the eigenvalues), which shows you that the two matrices must always have the same eigenvalues, no matter what dimension you take!

Now let's come back to the question of whether or not there are only four maximally entangled states (or even only four entangled states). A maximally entangled state is defined by the fact that its reduced density matrices are the unit matrix. This will let us prove the following theorem:

Theorem: Setting as above. Let $|\psi\rangle\in \mathcal{H}\otimes \mathcal{H}$ be a maximally entangled state, then there is a unitary matrix $U$ such that $|\psi\rangle =(1\otimes U)|\Omega\rangle$.

Proof: The idea is simple. We know by our observation that there exists a linear map $R$ such that $|\psi\rangle=(1\otimes U)|\Omega\rangle$ and we know how the reduced operators of $|\psi\rangle$ look like via the other observation. Putting this together, we obtain:

$$ \rho_1(|\psi\rangle\langle \psi|)=\sum_{i=1}^n \frac{1}{d}|R|f_i\rangle|^2 |i\rangle\langle i| \\ \rho_2(|\psi\rangle\langle \psi|)=\sum_{i=1}^n \frac{1}{d} R|f_i\rangle\langle f_i|R^* $$

(I recommend sitting down and writing out the details). In order for both matrices to be equal to $1/d$, you can see very quickly that the $R|f_i\rangle$ form an orthonormal basis, hence $R$ must be unitary. Conversely, for any unitary $R$, the two matrices are equal to $1/d$. End of Proof.

What about entangled states that are not maximally entangled? Well, as long as in the Schmidt decomposition, $k>1$, the state must necessarily be entangled, because it cannot be written as a product (this would be the case $k=1$).

Putting all this together, we can simply write down a state that is

a) a pure qubit x qubit state,

b) not maximally entangled but also not separable,

c) such that its reduced density matrices are not equal:

$$ |\psi_{sasha}\rangle = \sqrt{\frac{4}{9}} |0\rangle\otimes \left(\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)\right)+\frac{2}{3} |1\rangle\otimes \left(\frac{1}{\sqrt{2}}(|0\rangle-|1\rangle)\right)$$

Finally, what about mixed states? Can we get rid of the assumption that the eigenvalues must be the same?

The answer is yes. First let me note that the question of maximally entangled states is more intricate with mixed states. There is no obvious criterion. It's also a lot harder to work with, so I won't prove anything. So let's write down a state: Let $\rho_A,\rho_B \in \mathcal{B}(\mathcal{H})$ and $\rho_A\neq \rho_B$ such that they have also different eigenvalues. Then we can write down a state

$$ \rho_{sasha}=(1-p)|\psi_{sasha}\rangle \langle \psi_{sasha}|+p\rho_A\otimes \rho_B$$

As long as $p$ is small, this state should be entangled. Its reduced density matrices are given by $\rho_{1/2}=(1-p)\rho_{sasha,1/2}+p\rho_{A/B}$, and these two density matrices cannot have the same spectrum, because $\rho_{A/B}$ do not have the same spectrum.