  DAG models are statistical models satisfying a collection of conditional
independence relations encoded by the nonedges of a directed acyclic graph
(DAG) $\mathcal{G}$. Such models are used to model complex cause-effect systems
across a variety of research fields. From observational data alone, a DAG model
$\mathcal{G}$ is only recoverable up to Markov equivalence. Combinatorially,
two DAGs are Markov equivalent if and only if they have the same underlying
undirected graph (i.e. skeleton) and the same set of the induced subDAGs $i\to
j \leftarrow k$, known as immoralities. Hence it is of interest to study the
number and size of Markov equivalence classes (MECs). In a recent paper, the
authors introduced a pair of generating functions that enumerate the number of
MECs on a fixed skeleton by number of immoralities and by class size, and they
studied the complexity of computing these functions. In this paper, we lay the
foundation for studying these generating functions by analyzing their structure
for trees and other closely related graphs. We describe these polynomials for
some important families of graphs including paths, stars, cycles, spider
graphs, caterpillars, and complete binary trees. In doing so, we recover
important connections to independence polynomials, and extend some classical
identities that hold for Fibonacci numbers. We also provide tight lower and
upper bounds for the number and size of MECs on any tree. Finally, we use
computational methods to show that the number and distribution of high degree
nodes in a triangle-free graph dictates the number and size of MECs.
