Apart from telling that the hypergeometric functions (or series) are the solutions to the (essentially unique?) fuchsian equation on the Riemann sphere with 3 "regular singular points", the wikipedia article doesn't illuminate much about why this kind of special functions should form such a natural topic in mathematics (and in fact have been throughout 19th century). Simply:
What are hypergeometric series really, and why they should be (or have been in the past centuries) important/interesting?
In the 19th century, a lot of efforts were made in order to solve the general quintic equation $x^5+a_4x^4 +a_3x^3 +a_2x^2 +a_1x +a_0$ using special functions. It turns out that the roots of this equation are expressible in terms of hypergeometric series. To wit, one possibility is by first reducing the number of parameters, to the form $x^5-x-t=0$. Then a Lagrange inversion argument essentially gives a root $$ z=t {}_4 F_3(\frac15,\frac25,\frac35,\frac45,\frac12,\frac34,\frac54,\frac{5^5}{4^4}t^4)=t+t^5+10\frac{ t^9}{2!}+15\cdot 14 \frac{t^{13}}{3!}+\ldots $$
Hypergeometric series are solutions of a large class of differential equations. A series $\sum_{k} a_k t^k$ is hypergeometric if $Q_{k}=\frac{a_{k+1}}{a_k}$ is a rational function. Many familiar functions (trigonometric functions, exponential,logarithm,Hermite polynomials, Laguerre polynomials, etc) are hypergeometric.
Hypergeometric functions show up as solutions of many important ordinary differential equations. In particular in physics, for example in the study of the hydrogene atom (Laguerre polynomials) and in simple problems of classical mechanics (Hermite polynomials appear in the study of the harmonic oscillator).
Hypergeometric functions are also important in the study of elliptic elliptic curves where they can be used to compute the inverse of the $j$-invariant.
I guess you can read more about them in this wikipedia page or in these notes. Several examples of applications to number theory, physics and combinatorics can be read here .
One possible answer is that hypergeometric series were (and are) used to compute periods of elliptic integrals.
In modern terminology, take a smooth cubic $X \subset \mathbb{P}^2$ whose Weierstrass form is
$y^2w=x(x-w)(x-\lambda w), \quad \lambda \in \mathbb{P}^1-\{0, 1, \infty\}$.
Then $X$ is an elliptic curve, then it can be written as $\mathbb{C}/\Lambda$, where $\Lambda$ is a lattice. It turns out that the generators $\omega_1(\lambda)$, $\omega_2(\lambda)$ of $\Lambda$, i.e. the periods of the associated Weierstrass $\wp$-function
$\wp(z; \Lambda):=\frac{1}{z^2} + \sum_{l \in \Lambda-0} \big(\frac{1}{(z-l)^2}-\frac{1}{l^2} \big)$,
can be written in terms of the standard hypergeometric series $F$, namely
$\omega_1(\lambda)=i\pi F(\frac{1}{2},\frac{1}{2} , 1, 1-\lambda)$,
$\omega_2(\lambda)=i \pi F(\frac{1}{2}, \frac{1}{2}, 1, \lambda)$.
For further details see Chapter 1 of Kobliz's book "Introduction to elliptic curves and modular forms".
Hypergeometric functions arise as matrix coefficients of representations of Lie groups.
This is my formal answer, but I will also describe very informally why I believe this is an answer to your question. (By this I roughly mean: why this might have triggered interest of 19th century mathematicians who didn't have the language to be aware that this is what they were looking at. However: I do no not know a thing about the history of the subject so this is just a mathematical remark from a modern perspective.)
Interesting geometric spaces tend to have a rich group of symmetries (either because that is what makes them interesting (e.g. the circle), or because the symmetry is the only thing that gives us any grip on the object (e.g. spacetime)). Also studying a space through the functions on it has proved to be a powerful way to do mathematics. Hence representations of Lie groups on vector spaces of functions arise naturally.
Now it turns out that the matrix coefficient functions (of smaller representations) really give you a grip on these rather huge representations. The canonical example is when the group is compact and the interesting space on which the functions live is the group itself. In this case the Peter-Weyl theorem states that matrix coefficient functions form a basis. (Which at least makes it kind of credible that they are also useful for understanding functions on other compact homogeneous spaces.)
Now not everything is compact. (This is relevant here as Gauss' hypergeometric function appears as a matrix coefficient of an $SL(2, \mathbb{R})$-representation. (Incidentally this was the answer to prof. Koornwinder's question at my own thesis defense, see my comment to the original post above.)) However, also in this case matrix coefficient functions help you understand general representations on function spaces. For instance, in the infinite dimensional case (for non-compact groups irreducible representations need not be finite dimensional) it is not immediately clear (to say the least) that a given Lie algebra representation integrates to group representation. One way to make this work (given some extra conditions) is to first construct the matrix coefficients of the hypothesized group representation and then construct the actual representation from that (see for instance [1]). To me this always felt a bit like proving the existence of the Yeti from its footprint, but the difference is that it actually works.
Final remark: this is just the first part of JME's answer in disguise since in the setting of a Lie group acting on a space of functions on a homogeneous space (for that Lie-group), the differential equations JME is talking about come from the action of the Lie algebra.
[1] v.d. Ban: Induced representations and the Langlands classification [http://www.staff.science.uu.nl/~ban00101/manus/edinb.pdf]
The (general) hypergeometric equation has one more property which has not yet been mentionned: it has as (formal) solution at 0 exactly a series whose sequence of coefficients satisfies a first-order linear recurrence equation with polynomial coefficients. [One has to include all pFq here, including divergent ones].
The result above is a formal version of the action of the inverse Mellin transform on linear recurrence equations with polynomial coefficients.
In other words, structurally speaking, the hypergeometric equation is 'first order' (because its coefficient sequence is, not because its differential equation is), but it is the most general such.
The fact that this corresponds to a lot of known functions, as well as showing up in quite a few other places, is (often) a reflection of these structural properties.
You might like the wonderful book A=B:
