Fractal scaling--a power-law behavior of the number of boxes needed to tile a given network with respect to the lateral size of the box--is studied. We introduce a box-covering algorithm that is a modified version of the original algorithm introduced by Song [Nature (London) 433, 392 (2005)]; this algorithm enables easy implementation. Fractal networks are viewed as comprising a skeleton and shortcuts. The skeleton, embedded underneath the original network, is a special type of spanning tree based on the edge betweenness centrality; it provides a scaffold for the fractality of the network. When the skeleton is regarded as a branching tree, it exhibits a plateau in the mean branching number as a function of the distance from a root. For nonfractal networks, on the other hand, the mean branching number decays to zero without forming a plateau. Based on these observations, we construct a fractal network model by combining a random branching tree and local shortcuts. The scaffold branching tree can be either critical or supercritical, depending on the small worldness of a given network. For the network constructed from the critical (supercritical) branching tree, the average number of vertices within a given box grows with the lateral size of the box according to a power-law (an exponential) form in the cluster-growing method. The critical and supercritical skeletons are observed in protein interaction networks and the World Wide Web, respectively. The distribution of box masses, i.e., the number of vertices within each box, follows a power law Pm(M) approximately M(-eta). The exponent eta depends on the box lateral size l(B). For small values of l(B), eta is equal to the degree exponent gamma of a given scale-free network, whereas eta approaches the exponent tau=gamma/(gamma-1) as l(B) increases, which is the exponent of the cluster-size distribution of the random branching tree. Finally, we study the perimeter H(alpha) of a given box alpha, i.e., the number of edges connected to different boxes from a given box alpha as a function of the box mass M(B,alpha). It is obtained that the average perimeter over the boxes with box mass M(B) is likely to scale as <H(M(B))> approximately M(B), irrespective of the box size l(B).