Linear combo of normals is normal; how about other distributions? We know that univariate normal distributions are independent only if
their every linear combination is itself normal:
$$\tag{1}
  Z_i\sim N(\mu_i,\sigma_i^2)\,(\forall i)
  \implies\\ \biggl( \textrm{all $Z_i$'s mutually independent}
  \implies\Bigl(\alpha^T Z \sim N\bigl(\alpha^T\mu,\sum_i\alpha_i^2\sigma_i^2\bigr)(\forall \alpha\in\mathbb R^n)\Bigr)\biggr)
$$
I believe that the above holds for any $\mu\in\mathbb R^n$ and $\sigma\in\mathbb R_+^n$.
I am interested to know: does this hold for the normal distribution only?
Or are there other parametric distributions that have a similar property?
Suppose that $\{P_\theta\,:\,\theta\in\Theta\}$ is some family
of distributions, and that we have some $\theta_1,\theta_2,...\theta_n\in\Theta$.
When will the following be true?
$$\tag{2}
  X_i\sim P_{\theta_i}\,(\forall i)
  \implies \\ \biggl( \textrm{all $X_i$'s mutually independent} \implies
  \Bigl(\alpha^T X \sim P_{k(\alpha,\theta)}(\forall \alpha\in\mathbb R^n)\Bigr)\biggr)
$$
Here I have used $k:\mathbb R^n\times\Theta^n\to\Theta$ to denote some function that
generalizes the map
$(\mu,\sigma)\mapsto(\alpha^T\mu,\sum_i\alpha_i^2\sigma_i^2)$ from from the
normal case.
EDIT
After considering Ben's answer, I have changed the independence condition on the variables from "pairwise independent" to "mutually independent".
EDIT 2
Originally I was interested in a reverse implication
$$\tag{3}
  X_i\sim P_{\theta_i}\,(\forall i)
  \implies \\ \biggl( \textrm{all $X_i$'s mutually independent} \impliedby
  \Bigl(\alpha^T X \sim P_{k(\alpha,\theta)}(\forall \alpha\in\mathbb R^n)\Bigr)\biggr).
$$
See the comments on the accepted answer for an explanation of why this is not realistic.
 A: In your question you impose the condition of pairwise independence, but it is not possible to find the distribution of a linear combination of more than two random variables with this specification.  I will therefore assume that you are willing to impose the stronger condition of mutual independence, so that the distribution of the linear combination is fully determined.
On this basis, what you are looking for is the family of stable distributions, which can be obtained by considering the properties of characteristic functions (CFs) (i.e., Fourier transformations of the density functions).  If $X_1, ..., X_n$ are mutually independent random variables with characteristic functions denoted by $\varphi_k(t) = \mathbb{E}(\exp(it X_k))$, then any linear combination of these random variables has a CF that is a simple product of the CFs of the underlying random variables:
$$Y = \sum_{k=1}^n \alpha_k X_k \quad \quad \implies \quad \quad \varphi_Y(t) = \prod_{k=1}^n \varphi_k(\alpha_k t).$$
What you are looking for are classes of distributions where the form of the distribution is closed under linear combinations.  This property occurs for classes of distributions where the form of the characteristic function is closed under multiplication.  This occurs for the normal distribution as a special case, but it also occurs for a broader family of distributions called the Lévy alpha-stable distribution (see below).
Normal distribution: This has characteristic function $\varphi_k(t) = \exp(i \mu_k t - \tfrac{1}{2} \sigma_k^2 t^2)$.  We can see that this functional form is closed under multiplication:
$$\begin{equation} \begin{aligned}
\varphi_Y(t) = \prod_{k=1}^n \varphi_k(\alpha_k t) 
&= \prod_{k=1}^n \exp(i \alpha_k \mu_k t - \tfrac{1}{2} \alpha_k^2 \sigma_k^2 t^2) \\[8pt]
&= \exp \Big( \sum_{k=1}^n (i \alpha_k \mu_k t - \tfrac{1}{2} \alpha_k^2 \sigma_k^2 t^2) \Big) \\[8pt]
&= \exp \Big( i \mu_Y t - \tfrac{1}{2} \sigma_Y^2 t^2 \Big), \\[8pt]
\end{aligned} \end{equation}$$
where $\mu_Y = \sum_{k=1}^n \alpha_k \mu_k$ and $\sigma_Y^2 = \sum_{k=1}^n \alpha_k \sigma_k^2$.  From this characteristic function we see that the linear combination of random variables is also normal, with updated mean and variance that are corresponding linear combinations of the underlying means and variances.
Lévy alpha-stable distribution: This is a generalisation of the normal distribution, allowing for a generalised characteristic function that is still closed under multiplication.  It has characteristic function defined by:
$$\varphi_k(t) = \exp(i \mu_k t - \tfrac{1}{\gamma} |\sigma_k t|^{\gamma} (1-ib_k \cdot \text{sgn}(t) \Phi_{\gamma}(t))),$$
where $0 < \gamma \leqslant 2$ is a fixed parameter (common to each of the underlying distributions) and:
$$\Phi_\gamma(t) = \begin{cases} 
\tan (\tfrac{\pi \gamma}{2}) & & \text{for } \gamma \neq 1, \\[4pt] 
- \tfrac{2}{\pi} \ln |t| & & \text{for } \gamma = 1.
\end{cases}$$
(We also have the parameter restriction $|b_k| \leqslant 1$.)
This characteristic function is also closed under multiplication (when the parameter $\gamma$ is common to all the underlying random variables), and so it defines a broader class of distributions than the normal distribution.  The normal distribution occurs in the special case where $\gamma = 2$ and $b_k = 0$.  The corresponding density function for this distribution does not exist in closed form, but it is possible to derive many of its properties (which I will leave beyond the scope of the present answer).
